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SUMMARY 


This report traces the evolution in the use of state estimation for the analv- 
of aircraft flight data. A unifying mathematical framework for state estimation 
is reviewed, and several examples are presented that illustrate a general approach 

that'are^diff icult'To"' aCCUraCy data c °"sistency, and for estimating variables 
. t p . d ff ^ Xt fc measure. Recent applications associated with research- 
aircraft fiight tests and airline turbulence upsets are described. A computer 

LtPniwi T airCraft State estimation is discussed in some detail. This document is 

Aircraft ir * t ^' 3 manUal for the P r °g r am, called SMACK (SMoothing for 

h . 1 ^ ma X ° S ‘ e diversity of the applications described in the report 
emphasizes the potential advantages in using SMACK for flight-data analysis. 
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1 . INTRODUCTION 


Accurate determination of aircraft motions from noisy or incomplete measure- 
ments is an important problem in the analysis of flight-test experiments. The 
measurements often may contain significant errors which must be identified before 
the data are used in any performance or stability-and-control calculations. Fur 
thermore, direct measurements of certain important dynamic variables may be unreli- 
able or impractical to perform. A similar problem occurs in the analysis of air- 
craft accidents, where the actual motions may have to be determined from a very 
limited data set. These problems are being solved by the analytical method known as 
state estimation. This report presents an algorithmic approach for aircraft state 
estimation, demonstrates its application for solving several example problems, and 
describes the computer program used to obtain the solutions. 


The first application of state estimation to postflight data analysis can 
probably be attributed to the pioneering work of Otto Gerlach m the 1960s at the 
Delft Technological University, The Netherlands. This early contribution (refs, 
and 2) called "flightpath reconstruction," was primarily concerned with the accu- 
rate determination of angle of attack, pitch angle, and vehicle velocity during 
dynamic maneuvers. These "states" were obtained by integrating functions of mea 
surements from the pitch-rate gyroscope and normal and longitudinal accelero- 
meters Initial conditions and bias terms were determined from airspeed and alti- 
tude measurements at steady-state e-d points of the maneuver. The resulting 
"smoothed" time-histories were then used as a basis for subsequent parameter identi- 
fication studies. 

Application of the state-estimation method to aircraft problems is possible 
because the forces and resulting motions of an aircraft along a flightpath are 
related by well-known equations of motion. The equations may be used to produce 
estimates of force and motion variables that are compared with corresponding mea- 
surement time-histories in an iterative procedure until a suitable "match" is 
obtained. As Gerlach has pointed out, the technique of state estimation provides a 
check on instrument accuracy and data consistency as well as estimates of unmeasured 
or poorly measured variables. 

These items have been the primary objectives in most of the studies that fol- 
lowed the initial work of Gerlach. His students later improved and formalized the 
techniques that Gerlach had developed (refs. 3 and 4). In this country, early 
advocates of the use of state estimation for flightpath reconstruction were Wingrove 
(refs 5 and 6) at NASA Ames, Eulri.'h and Weingarten (ref. 7) at Calspan, an 
Molusis (ref. 8) at Sikorsky Aircraft. Over the past few years, the work in this 
field has been evolving consistently toward the development of more sophisticated 
algorithms, the use of more complete kinematic models, and the treatment of more 
difficult applications. 
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Most recently-developed algorithms (refs. 9-18) utilize a version of 
extended Kalman filter. Although good results have been reported in offline flight 
data-processxng applications, such algorithms are not optimum in their of 

aircra^statfestimar ^ ^ aSU ^ emenfc record - The algorithm advocated here for 

Til ? c : a ZlTT t ° f a nonlinear - fi * ad - 

estimates until a minimum squared-error measure is achieved ” Pr °' eCl state 

a nominal trajectory and convergence is quadratic uT h’ blnearlzatlo n is about 

sweep" algorithm of HcReynolds and Bryson (re? n, , ,!? 7 

optimal control problem This alfforit-hm h ' u finally devised to solve an 

aircraft motions along a flight !'f ' had n0t been a PP lletl to determine 

mation program des or bed X rZ T of “» «ti- 

Kinematics, called SMACK, ^e^Le^d aTE’ t 

0? a al??? a ?? a^i?;??; 1 the Natl0nal T ™ s P°n ta “°" safety Board in its investigation 

effecute^ oTthTslfptogX* tVZttuT ‘"1°™??" 

ings, includes discussions of the underlying mathematics" ?“ al * obaptar head - 
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forward smoother is shown to have^ °° nS1StS ° f a backward information filter and 
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2. STATE-ESTIMATION ALGORITHMS 


In this report, the tehm "state estimation" refers to offline processing of a 
set of independent measurements from a given physical system. Each record in the 
set consists of a time history that covers the same fixed interval, b!t does noT 
necessarily share a common sampling rate with other records in the’set If the 

m^icaUy an^he g meas ^1° ° f the s y stem are understood mathe- 

y, all the measurement records may be processed together in an "optimal" 

.ay. m this case, the objective is to determine a set of initial conditions and 

cing functions that .ill cause the output of a mathematical model to "match" the 

measurement time histories, usually in a least squared-error sense The “ „ r 

postflight state estimation is known as a fixed-interval smoothing problem WUh 

s solution, the analyst can determine that the measurements are consistent and in 

errorn^caWactors'!^ 83 ^ tables, as well as instrument bias 

The state-estimation process solves a state model 


( 2 . 1 ) 


x(i + 1) = f[x(i) ,w(i) ] ; x(0) = x Q 

such that y(i + 1) in the measurement model 

z(i + 1) = y(i + D + v(i -*■ 1) ; y ( i + 1) = h[x(i + i)] {2 .2) 

acceptably matches the data record over a time interval (i = 0,1,- N - 1) usuallv 
n a least-squared error or minimum variance sense. In equation (2 1) x(i) is an 
State -f- /„d w(i) is an NW-element forcing-functiot lector i“ “ 
meets a ’ * V( * ” are vectors representing the measure- 

r and N r:;rr measurement — , y<1 . „ ls t j syste :zz t 

variables ° PU Ve0t ° r iS « enerall y a nonlinear function of the state 

finftl Hfrr ra ^ C appllcatlons ' the state and measurement models together represent 

bod? As desr re h C H apP I PXl ” atl0 " f ° r the ^-“egree-of-freedom dyniics o 
body. As described in the next chapter, the models are used to generate time hi! 

ries which are likely to be found in a flight-test measurement set These include 
onboard variables such as Euler angles, angular rates, and linear accele^atLns as 

er“rs S or scal'e^rt 1315183 SU ° h “ Slant h 3 " 6 "’ bearIng ’ a " d Nation. Any biis 
to the Z factors associated with the state or measurement models are appended 

to the state vector and treated as constant but unknown parameters. PP 

the . Thas ° luti ° n of the fixed-interval smoothing problem consists of determining 
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(2.3) 


J = (-0 - V tp c ’<*0 


x )/2 

o 


N-1 


I 


i=0 


{[z(i + 1) - y(i + 1)] t R" 1 [z(i 


1) 


y ( i + 1)] + w fc (i)Q 1 w(i)}/2 


estimate of x n , and P Qt Q, and R are 


is an a priori estimate ui a q , auu i 0 , 

°Sage and Melsa (ref. 21) discuss a maximum-likelihood interpre- 


In equation (2.3), x 

tation of the performance measure in which P Q is the error-covariance matrix for 
the a priori estimate, and Q and R are error-covariance matrices for the input and 
output sequences, respectively (assumed to be stationary). Note that the first term 
of equation (2.3) serves as a "penalty" function and tends to bias the estimate of 
x Q toward its a priori value. 

The fixed-interval smoothing problem is solved using a method of successive 
approximations based on expansion of the performance measure (eq. (2.3)) to the 
second order and of the dynamic constraint (eq. (2.1)) to the first order Suppose 
we choose x , {w(i)J and obtain a nominal trajectory by solving equation (2.1). It 
is unlikely that our solution minimizes J, but we shall try to determine a neigh- 
boring solution that yields a smaller value. To do this, we first express a varia- 
tion in the performance measure in terms of variations 6x Q and l«w(i )}: 


N-1 


6J = (x Q - x 0 ) tp 0 1|5x 0 + (1/2)6x 


p-v-Y 

° ° ° Lu 


{w t (i)Q _1 6w(i) + ( 1/2)6w"(i)Q 


t, . , A -1 


i=0 


6 w(i) - [z(i + 1 ) - y(i ♦ 1)]V 1 h x 6x(i + 1 ) + ( 1 / 2 ) 6 x"(i ♦ DhjR'V* 


t/. 


X (1 +1)} 


(2.4) 


Next, we assume that deviation 

6x( i + 1 ) = 

In equations (2.4) and (2.5), 
defined as 

f x = 3 f[x(i),w(i)]/ax(i) ; f w 


from the nominal trajectory will be governed by 
f 6x( i) + f u 6w(i) , «*(0) = 6 X 0 (2>5) 

X w 

the Jacobian matrices (partial derivatives) are 

= 3f[x(i),w(i)]/3w(i) ; h x = 3h[x(i + 1)]/3x(i + 1) 

( 2 . 6 ) 


and are to be evaluated along the nominal trajectory. 

Our objective now is to specify 6* 0 and UwU» such that U has .the most 
negative value possible, subject to the dynamic constraint of equation (2.5). We 
solve this "accessory minimization” problem In the usual uay (ref. 22) by adjoining 
the constraint to equation (2.4) using a Lagrange multiplier. Hence, 
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(2.7) 


6 J 


N-1 


6 J + 



i=0 


x fc ( i + 


1) 


f 6x(i) + f 6w(i) - 6x( i + 1)] 

A W 


The necessary conditions for minimizing J lead to a linear, two-point boundary- 
value problem (LTPBVP) given by equation (2.5) and 

*(i) = ^(*(1 + 1) + h^R ^h x 6x(i + 1) - h^R~^[z(i + 1) - y(i + 1)]} 

X(°) = -p; 1 [x 0 - x 0 ) + 6x o ] ; X(N) = 0 (2.8) 

6w( i ) = -w(i) - QfV fc x(i) (2 9) 

This LTPBVP has an exact solution (ref. 22). Hence, it is possible to determine 
6x q and { 6w ( i ) } , recompute the nominal trajectory with 

X o * x 0 + 6x 0 J w (i) * w(i) + 6w( i) 

and evaluate the performance measure, iterating until J is minimized. The change 
in J that should be realized at any iteration is found by substituting equa- 
tions (2.8) and (2.9) into equation (2.7). 


6 J 


N-1 

^ ^ [<Sw t (i)Q ^wfi) + 6x t (i + 1 )h^R _1 h x 6x( i + 1)]/2 (2.10) 
i=0 


Two equivalent sweep solutions of the LTPBVP are given here. The first is 
derived by introducing a vector 6x(i) and matrix P(i) and letting 

6x(i) = 6x(i) - P(i)x(i) (2.11) 

Notice that the boundary conditions of equation (2.8) require that 

6x(°) = x Q - x Q ; P(0) = P 0 ; 6x(N) = «5(N) 

Straightforward algebraic manipulation yields the algorithm outlined in 
table 2.1(a), which is essentially the procedure proposed by Cox in 1965 
(ref. 23). It consists of a forward covariance filter and a backward smoother a 
form that invites comparison with the extended Kalman filter often employed for 
nonlinear state and parameter estimation (ref. 24). We observe that, for a class of 
systems with a state model that is linear in its forcing function, as 

f[x(i),w(i)] = g[x(i) ] + Gw( i) 
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TABLE 2.1.- EQUIVALENT ALGORITHMS FOR SOLUTION OF A NONLINEAR SMOOTHING PROBLEM 

With x and {w(i)} obtained from the preceding iteration (or an initial guess), 
compute a nominal trajectory using equations (2.1) and (2.2) and evaluate the 
performance measure using equation (2.3)- Now perform (a) or (b) as follows. 


(a) Forward filter/backward smoother (b) Backward filter/forward smoother 


Forward filter with a time update 

6x( i) = f 6x(i - 1) - f w ( i - 1) 

X w 

HU) = f/U - i)f * 

5x(0) = x q - x q ; P(0) = P Q 

and a measurement update 
6x( i ) = 5x( i ) + K(i)e(i) 
p(i) = [I - K(i)h x ]M(i) 
where 

e(i) = [z( i) - y(i) ] - h x 6x(i) 

K( i) = M(i)h x R _1 
R = R + h M(i)h fc 

X X 

Backward smoother 

S(i) = [I - K(i)h ] fc [x( i) - hV 1 e(i)] 

X A 

X ( i - 1) = f x B(i) ; UN) = 0 
and 

6w( i - 1) = -w(i - 1) - Qf^B(i) 


Backward filter with a measurement 
update 

B(i) = a( i ) - h x R _1 [z(i) - y(i)] 

S(i) = M(i) + 
a(N) = 0 ; M(N) = 0 
and a time update 
a(i - 1) = f fc [B(i) - S(i)f d(i)] 

X W 

M(i - 1) = f fc [I - f L(i)] fc S(i)f 

X w A 

where 

d(i) = Q[f fc B(i) + Q _1 w( i - 1)] 
w 

L(i) = Qf fc S( i ) 

W 

Q = [Q- 1 * r‘s(i)fj-' 

Forward smoother 
Sx(i + 1) = f 6x(i) + f w «w(i) 

6x(0) = 6 x q 
where 

6x = -[P' 1 + M(O)] -1 
0 o 

- mo) * p;'(x 0 - 

6w(i) = -d(i + 1) - L(i + l)f x «x(i) 


o o 


Update x and (w(i)}. Loop and iterate until 6x Q and (6w(i)} are "sufficiently" 
small and° the performance measure is minimized. 
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the forward covariance filter of table 2.1(a) is identical to an extended Kalman 
filter linearized about a (prior) nominal solution. The usual linearization, how- 
ever, is about a current solution. In at least one case (ref. 25), the extended 
Kalman filter has been coupled with a backward smoother. Such a procedure requires 
no starting solution but does not iterate to minimize a performance measure, and so 
provides only an approximate solution of the nonlinear smoothing problem. 

A second, more useful sweep solution of the LTPBVP is obtained by introducing a 
vector a(i) and matrix M(i) and letting 

X(i) = a(i) + M(i)<$x(i) (2.12) 

In this case the boundary conditions of equation (2.8) require that ct(N) = 0 
M(N) = 0, and 

6x q = -[P^ + M(0)]' 1 [P^(x o - x q ) + a(0) ] (2.13) 

The resulting algorithm is outlined in table 2.1(b) and consists of a backward 
information filter and a forward smoother (for offline application). It can be 
shown that this algorithm is equivalent to the "modified" Newton-Raphson method if 
there are no unknown forcing functions (ref. 26). Notice that the sequences (d(i)}, 
(L(i)}, computed during the filter pass, are utilized during the smoothing pass. 

Here the (temporary) storage requirement depends on the dimensions of w(i) and x(i) 
and, of course, on the length of the data record. This formulation, which has been 
implemented with SMACK, has the following advantages over the algorithm of 
table 2. 1(a) : 

1. The a priori covariance is easily specified by setting P' 1 = 0 in equa- 
tion (2.13), which is often a good choice in practice. This is equivalent to "no 
confidence" in an a priori estimate. 

2. Constant elements (bias-error and scale-factor parameters) of the state 
vector are naturally decoupled from the dynamic states. This feature reduces the 
computational burden. 

In applying either algorithm of table 2.1, the analyst must be careful in 
choosing starting values for x Q and { w( i ) } , and in selecting the weighting 
matrices P Q , Q, and R. The convergence properties of the algorithm are influenced 
directly by the nominal solution generated by the initial choice of x and any 
unmeasured forcing functions. Suitable starting values can be obtained by solving 
the finite-difference approximation of the state model for the forcing-function 
sequences by using filtered versions of the measurement records to construct the 
state estimates. On the other hand, the nature of the solution, once convergence is 
obtained, depends to a considerable degree on the choice of the weighting matrices 
p o» Q* a £ d R * Th e effect of P Q is to bias the estimate x Q toward the a priori 
value x q ; it can usually be ignored when applying the algorithm of table 2.1(b) 
Reasonable values for the elements of Q and R may be determined as follows: 
filter each measurement record until the residual sequence appears sufficiently 
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"white," and use its variance as the appropriate diagonal element of R; then con- 
struct an estimate of each forcing function, and use the mean-square value of the 
starting sequence (w(i)} as the appropriate diagonal element of Q. 

It should be observed in passing that, heuristically , smoothing is a process of 
zero phase-shift filtering in which bandwidth increases as the scale of Q 
increases and the scale of R decreases. One expects that forcing-function and 
residual variances will agree with the corresponding elements of Q and R used in 
obtaining the solution. That solution, however, is not unique. Scaling the 
elements of Q up and the elements of R down by the same factor (bandwidth 
increase) will result in a solution having closer fits to the data, but with 
"noisier" forcing-function estimates. This situation emphasizes the need for the 
analyst to carefully consider the engineering aspects of the problem. 

Lest the potential SMACK user despair over possible pitfalls to be encountered 
in choosing starting values and selecting weighting matrices, he or she should be 
assured that the program has been designed to require little user intervention. A 
subroutine included with SMACK provides the set of initial conditions and forcing 
functions needed to generate a starting trajectory. This routine also calculates 
sets of diagonal element values for the Q and R weighting matrices, using the 
procedure suggested earlier. An outline and block diagram for the starting 
subroutine are given in appendix A. 


The Linear Case 

Finally, for the sake of completeness, a formulation of the algorithms of 
table 2.1 for a linear system is considered that is useful in offline digital fil- 
tering applications. For the linear case, the state and measurement models of 
equations (2.1) and (2.2) become 

x(i) = Fx( i - 1) + Gw( i - 1) , x(0) = x Q , z(i) = Hx(i) + v(i) (2.14) 

where F is an NX*NX matrix, G an NX*NW matrix, and H is an NV*NX matrix. Using 
the notation expressed in equation (2.6), we obtain for the Jacobians 

f x = F ; f w = G ; h x = H (2.15) 

Now consider a simple change of variable 

6x(i) = x(i) - x p (i) (2.16) 

where x (i) is any (nominal) solution of equation (2.14) and x(i) is to be the 
solution of the linear fixed-interval smoothing problem. If equations (2.15) 
and (2.16) are used in the forward-filter, backward-smoother algorithm of table 2.1, 
along with 


6x( i) = x(i) - x n ( i ) ; 6x(i) = x(i) - x n (i) 
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the algorithm of table 2.2(a) results. In similar fashion, if equations (2 15) 

and (2.16) are used in the backward-filter, forward-smoother algorithm of table 2 1 
along with * ’ 


a(i) <- a(i) + M(i)x n (i) ; e(i) <- 8(i) + S(i)x n (i) 
the algorithm of table 2.2(b) results. 

TABLE 2.2.- EQUIVALENT ALGORITHMS FOR SOLUTION OF A LINEAR SMOOTHING PROBLEM 
(a) Forward filter/backward smoother (b) Backward filter/forward smoother 


Forward filter with a time update 

x(i) = Fx( i - 1) 

M(i) = FP( i - l)F fc + GQG t 

x(0) = x ; P(0) = P 

o o 

and a measurement update 
x(i) = x(i) + K(i)e(i) 

P(i) = [I - K(i)H]M(i) 
where 

e(i) = z(i) - Hx( i ) 

K(i) = M(i)H t R" 1 
R = R + HM(i)H t 
Backward smoother 

8(i) = [I - K(i)H] fc [x(i) - H fc R _1 e(i)] 
X( i - 1) = F t 6( i) ; x(N) = 0 
and 

w(i - 1) = QG t 6(l) 

X o = ' P o X(0) 


Backward filter with a measurement 
update 

8(i) = a(i) - H t R" 1 z(i) 

S(i) = M(i) + H fc R -1 H 
a(N) = 0 ; M(N) = 0 
and a time update 
a(i - 1 ) = F t [s(i) - S ( i ) Gd ( i ) ] 

M(i - 1) = F fc [ I - GL( i ) ] t S( i )F 
where 

d(i) = QG fc 8(i) 

L(i) = QG fc S( i ) 

Q = [Q' 1 + G t S(i)G]" 1 
Forward smoother 
x(i + 1) = Fx( i ) + Gw( i ) 
x(0) = x Q 
where 

x 0 = ~ [P o ] + M (0)] _1 [“(0) - P -1 x ] 
u ° o o 

w(i) = -d(i + 1) - L(i + l)Fx(i) 


Each of these algorithms converges in one step: no starting solution is 

needed However, the considerations concerning the choice of weighting matrices Q 
and R, and the a priori estimate x Q , P Q are the same as for the nonlinear case. 
Hence the backward-filter, forward-smoother algorithm also has an advantage in 
solving any linear fixed-interval smoothing problem. That algorithm has been used 
in the realization of a low-pass filter employed in the SMACK starting procedure. 
The filter is described in appendix E. 
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3. STATE AND MEASUREMENT MODELS 


In this chapter the mathematical models utilized by SMACK to obtain flight 
rajectories are defined and discussed. Aircraft motions are assumed to be governed 
by a six-degree -of -freedom kinematic model, referred to a flat, nonrotating Earth 
The usual choice of state variables leads to a formulation in which both state and 
measurement models are nonlinear (ref. 25). The solution algorithm outlined in the 

matrices T 7 "and' BValu ^\ b ° th the state and measurement model Jacobian 
. f x’ f w’ and h x ln ec l* along the trajectory. However, if the state 

ear 1 ? ^ <*“"» *• to obtain a sLte model (hat L u„ 

ear. All nonlineanties then appear in the measurement model. The advantage of 

s'tate SMACK iS com P u ^ at i°nal : the Jacobian matrix X the 

the Pff? 13 TL al ° ng any trajectory, a feature that significantly improves 
the efficiency of the solution algorithm (ref. 27). 

. „. F ° r reallzafc ion of a linear state model, we start with the vehicle attitude 
defined by the Euler angles <♦, e, *), and the vehicle position, defined by the 

consis't U of a tL 0 O d° a t teS <X ; y ’- h) . as »!H*t variables. Other (bate variables 
consist of time derivatives U, 0 , *), (*, 0 , 4 ,), and (i, y, h), ( x v £) When 

velocuv states 1 ! C ° nSldered ' the state "°“el is augmented with wind- 

evide„t y in dLg^'sho^ 

b rSf 1 S 28 )° f Note 8 ^at r th bl, f dleS ’'' r Stru0ture knoun as a Brunovsky canonical form 
a w that the forcin S Unctions for this system are (d , d d u ) (d 

mnHei n • i h r g *’ §y ’ gh ^ ’ In S ° me s i !r uat ions ifc may be preferable X to use a simpler 
model with forcing functions (x, y, h) instead of (d , d , d K ) and/or T 7 

instead of (d. . d d ) x * h '’ ana/or (<t>, 0 , 41 ) 

m * n 7 * 

. , The measurements available in an aircraft state-estimation problem often 

linen -t- 6 ’ vehlole attICude . velocities, and accelerations. All non- 

linearities associated with aircraft kinematics appear in the measurement model 

shown on the right side of figure 3.1. For example, the blocks labeled with ^ £ 
represent the transformation from Earth-surface axes to vehicle-fixed body ax“ 

atio ns" a ^ a" f l6S “ Ib ? r6SPe0t t0 the alr ” aSS <“• *• “> -d Ihe tdy cel- 

erations (a x , a y , a z ) are calculated from 3 


= L 


x - w. 


y - w. 


n - Wu 


= L 


h + g 


(3.1) 


respectively, where the transformation is defined by the direction- 


cosine matrix 


13 












cos 9 cosi> 

cose simp 

sine 

sin<j> sine cosij) 

sinip sine simp 

-sin<|> cose 

- cos<() sim(j 

+ COS(p COSlp 


cost)) sine cosip 

cosip sine simp 

-COSt)) cose 

+ sin<)> simp 

- simp cosip 



(3.2) 


Radar measurement variables R (slant range), B (bearing angle), and E (eleva- 
tion angle), which are given by 


R = [(x - x p ) 2 + (y - y r ) 2 + (h - h p ) 2 ] 1/2 

B = tan" ^ [ (y - y r )/(x - x r )J ; E = sin -1 [(h - h r )/R] (3.3) 


where (x r , y r , h r ) is the tracking antenna location, are calculated in the block 
labeled with a 72 in figure 3.1. Other inertial data might be supplied by an 
onboard inertial navigation system (IMS). Although not shown in figure 3.1, the 
measurement model does include the INS velocity variables V (groundspeed) and <j> 
(groundtrack) , which are g g 

V g = + y 2 ) 1/2 > ♦g = tan _ 1 (y/x) (3.4) 

In order to fit air-data records of true airspeed and the flow angles, the 
measurement model can provide estimates of the aerodynamic variables V (airspeed), 
a (angle of attack), and 8 (sideslip angle), from the relations (ref. 29 ) 

V = (u 2 + v 2 + w 2 ) 1/2 ; a = tan _ 1 (w/u) ; 6 = tan _ 1 (v/u) (3.5) 


These are calculated in the block labeled .A. Notice that the variable 8 models 
the vane flow angle, which is not quite the same as the usual sideslip angle 
(ref. 30). When air-data measurements are included in the SMACK estimation proce- 
dure, the winds along the flightpath can also be estimated. The wind variables in 
the measurement model, W xy and W hd (horizontal magnitude and heading), and V 
(vertical magnitude) are calculated in the block labelled W from the relations 


u / 2 2.1/2 

W = (w + w ) 
xy x y’ 


W, 


hd 


_ 1 

= tan (-w /-w ) 

J X 


wd 


= w L 


(3.6) 


Blocks labelled H and Q represent the nonlinear relations that express the 
body angular velocities (p, q, r) and angular accelerations (a^, a^, a ) in terms of 
the state variables. The angular velocities are calculated from (ref. 30) 


p = 4 » - 4 » sine , q = e cos<t> + tp simp cose , r = -e simp + ip cos<p cose ( 3 . 7 ) 


and the angular accelerations are calculated from 
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a £ = <p - ip sin9 - 9ip cos9 

a = ill sin* cos9 + 9 cos4> - 9ip sin* sin9 + <pr 

m 

a n = 4> cos4> cos9 - 9 sin4> - 0ii> cos<t> sin9 - <t>q (3.8) 

If required, position corrections for location of air-data or accelerometer 
instruments can be made in the SMACK measurement model. Body-velocity corrections 
for the air-data system are given by (ref. 3D 

Au = qz i - ry i , Av = rx t - pz L , Aw = py t - qx L (3-9) 

where (x-, y-, z^) are the body coordinates of the instrument position with respect 
to the aircraft center of gravity. Airspeed (pitot-static) and aerodynamic angle 
(vane) corrections are treated separately. Corrections to the body accelerations 
are (ref. 3D 

Aa x = (pq - a n )y i + (pr + a m )z i - (q 2 + r 2 )x i 

Aa y = (pq + a n )x i + (qr - a^^ - (p 2 + r 2 )y t 

o o 

Aa z = (pr - a m )x i + (qr + a^^ - (p* + q )z t (3.10) 

Note that equations (3.9) and (3.10) refer the corrected variables to the instrument 
location. Furthermore, equation (3.10) utilizes angular acceleration estimates (a^, 
a , a ), which from equation (3.8) are seen to be functions of (<p, 9, ip). The 
algorithm of chapter 2 as implemented in SMACK requires all estimates formed in the 
measurement model to be functions of state variables. In this case, then, the user 
must specify the forcing functions to be (d^, d m , 8 n ). 

A general rule in the application of SMACK is that if there are any elements of 
sets (a , a , a ) or (a a , a m , a n ) to be estimated, then (d x , d y , d h ) or (d^, d m , d n ) 
must be specified as forcing functions. Otherwise, the forcing functions may be 
chosen from the sets (<t>, 9, ip) or (x, y, h). Care should be taken not to mix 

elements of (d^, d m , d n ) with elements of (<P, 0, 4») as forcing functions, or, 

similarly, elements of (d x , d , d^) with elements of (x, y, h). It is true, of 
course, that the more integrators there are between input and output of the state 
model, the more "smooth" the output will be. 

In certain situations, such as performing a preliminary data-consistency check, 
a user may wish to employ measurements of (a x , a y , a z ) and/or (p, q, r) to generate 
forcing functions. For example, accelerometer measurements can be used to produce 
the earth-frame accelerations required for the SMACK state model by solving equa- 
tion (3.D to obtain 
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In similar fashion, it can be seen that rate-gyro measurements will provide the 
required Euler-angle derivatives by solving equation (3.7) to obtain 


i = (q sin* + r cos*)/cos0 , 0 = q cos* - r sin* , * = p + * sin© (3.12) 

Clearly, if either equation (3.11) or (3.12) is employed, the state model will no 
longer be linear. The tradeoff is, of course, that time-histories for (x, y, h) 
and/or (*, 0, *) need not be estimated. A separate consideration here is that the 
use of noise-contaminated measurements to generate forcing functions is likely to 
bias the solution in an unpredictable way. 
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4. APPLICATION EXAMPLES 


State estimation as a means of checking instrument accuracy and data consis- 
tency is now used by many flight-test groups (refs. 9-18). Once a consistent, 
smoothed set of time histories is obtained from the data, other analyses, such as 
identification of stability and control derivatives, are readily performed. In 
fact, relatively simple routines may be used for identification tasks, allowing the 
analyst freedom to develop a proper aerodynamic model. Since the data-consistency 
application is extensively treated in chapter 6, and a flight-test example of a 
complete SMACK solution is given in appendix C, it will not be discussed here. 
Instead, some of the more recent applications of aircraft state estimation in 
obtaining estimates of unmeasured or poorly measured variables will be addressed. 

In this chapter four examples, based on recently-reported applications of 
aircraft state estimation, and one example of an application not previously 
reported, are discussed (ref. 32). The applications, quite diverse in terms of the 
available measurements and desired estimates, illustrate the wide range of problems 
that can be treated in a unified way by using SMACK. Data for each example were 
taken from a simulated trajectory consisting of a rising, coordinated, 180° turn m 
the presence of wind. The trajectory is generated by a SMACK subroutine for user 
testing of a problem coding list. Small amounts of random noise, usually \% or 
less, were added to each measured variable, and all measurements were recorded once 
per second. A summary of the available measurements and variables to be estimated 
for each example is given in table 4.1. 


Example 1 

For aircraft accident analysis, state estimation can be effectively used to 
combine data from several sources (e.g., radar site and flight recorder) to deter- 
mine motions along a trajectory (refs. 33 and 34). In addition, the winds along a 
flight trajectory can often be estimated. Wind estimation has been used m the 
analysis of recent airline turbulence upsets, and is the subject of the first 
example This application is covered in some detail in chapter 7. Parks et al . 

(ref. 35) describe the estimation of winds by using data from a DC-10 encounter with 
severe high-altitude turbulence. The wind estimates from that analysis led Parks to 
hypothesize the presence of a classical "cat's-eyes" vortex phenomenon in the jet- 
stream shear layer at the time of the encounter. 

Data from a digital flight recorder like the one carried by a DC-10 includes 
accelerations, Euler angles, altitude, and airspeed, sampled at intervals of 
0 25-4 0 sec Sufficient additional information is available to approximate the 
records of angles of attack and sideslip. The addition of ground-based air traffic 
control (ATC) radar provides a number of measurements approaching that available 
from flight test. To obtain the desired wind estimates, Parks first transformed the 
accelerations into an Earth frame, then integrated them to obtain aircraft velocity 
with respect to the Earth. A consistent set of initial conditions and accelerometer 
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TABLE 4.1.- LIST OF VARIABLES FOR STATE-ESTIMATION EXAMPLES 


Example 


Variable 

1 

2 

3 

4 

5 

Linear Acc. 






[ a x» a y ’ a z ^ 
Angular Vel. 

Measured 

Measured 


Measured 

Measured 

(p, q, r) 
Position 


Measured 



Measured 

(R, B, h) 
Winds 

Measured 


Measured 

Measured 

Measured 

( W X y> W hd» V wd^ 
Angles 

Estimated 

Measured 

Measured 

Measured 

Measured 

( 4> , 0, ip ) 
Air Data 

Measured 

Estimated 

Measured 

Measured 

Estimated 

6) 

Measured 

Measured 

Estimated 

Estimated 

Estimated 


^ ° btained b y etching calculated-position time histories with 
radar and barometric altitude records. The wind components were then found as the 

the aircraft veiocities with respect to E - th - a - f n he 


The first example illustrates a wind-estimation application and uses the mea 

upsensertable^^n 6 'i ^ anal * sis of the DC ~10 turbulence 

p t (see table 4.1). In the analysis of this problem by SMACK, all elements of 
the forcing-function vector (», e, *), (g x , g g ), and ( d d ’ dJ ar f‘ 

estimated. All of the measurement time histories (a a a ) y /u h Q \ . . 

B, h) ace fitted in the ieast-sduaned-ennor vroZlli’. fhe ^sulung’wiM “titles 
are shown in figure 4.1, along with the "true" winds for comparison. The close 
agreement of the horizontal wind records indicated there and in table 4.2 L prob- 
ably better than could be expected in practice, since ATC enroute radar data are 
recorded only about once every 10 sec. 


Example 2 

Other applications of state estimation that are becoming increasingly imnnrt-ant- 
are associated with the testing of high-performance aircraft In large angle-o? 

c maneuvers and spin tests, for example, measurements of Euler angles air 
speed and aerodynamic angles (e.g., angles of attack and sideslip) may contain' 
sig^ficant errors. In a recent paper, Taylor (ref. 36) discussed the estimation of 
u er-ang e ime histories and air-data instrument bias errors and scale factors for 

erometei" 8 T' " The measurement set for this application consisted of accel- 

known Wi[h th 8yr °’ and / lr ‘ data retirements. The winds were assumed to be 

Tailor "? tLd" aCCelerations and angular velocities as forcing functions, 

y tted the air-variable measurement records, using a squared-error 
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ESTIMATION OF WINDS 

ooooo TRUE SOLUTION 


ESTIMATE 




Figure 4.1.- Winds for example i. (a) Horiz. magnitude, (b) horiz. heading, 
u (c) vertical wind. 

criterion and a Newton-Raphson algorithm to determine the desired estimates of 
biases scale factors, and Euler-angle time histories. To avoid possible singu- 
larities in angle calculations, Taylor utilized the differential equations relating 
the angular velocities and direction cosines. 
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TABLE 4.2.- RESULTS FOR STATE-ESTIMATION EXAMPLES 


Variable 

i 


2 

3 


4 


5 


Winds 

Mean 

S.D. 

Mean S.D. 

Mean 

S.D. 

Mean 

S.D. 

Mean 

S.D. 

W xv' kt 
W hd’ de g 
V wd’ m/s 

0.03 

-0.84 

-0.01 

0. 12 
1.45 

0.03 








Angles 
<t>, deg 
0, deg 
i>, deg 



-0.05 0.10 

0.03 0.04 

-0.75 0.88 





0.00 

-0.04 

0.24 

0.07 

0.09 

0.29 

Air Data 
V, kt 

а, deg 

б, deg 




0.08 

0.02 

0.04 

0.44 

0.11 

0.16 

0.02 

0.01 

-0.02 

0.05 

0.05 

0.06 

-0.05 

0.00 

-0.22 

0. 16 
0.07 
0.28 

Note: Mean 
estimated 

and s.D. ( stanc 
time histories. 

lard deviation) 

refer 

to error between true and 

•- 


ample 2 illustrates the application of state estimation for determining Euler 
anges using the measurement set of Taylor as summarized in table 4 1 i n the 
analysis by SMACK, the inertial wind components (w , w W J we re obtain^ f „ 
measured winds, and used in equation (3.1). Estima^d lleZelZr the forcing func! 
were (<t>, 9, <4) and (d x , dy, d h ). The measurements (a Y , a , a i (n n r ) and 

Lr dlt 6 W6re fitt ®^' with bias-error and scale-factor estimates obtained^or the 
-data records. The Euler-angle estimates are shown in figure 4 2 plotted with 
the corresponding true values. Estimation errors are given in tabie\ 2 lrt 

e noted that the pitch-angle excursion is not large along the simulated ‘ trljec 
tory; For extreme maneuvers, in which the pitch angle may approach 90° it is not 

possible to avoid singularities using the linear (coordinate-transformed) staL 
model used in the SMACK program. ransrormecu state 


Examples 3 and 4 

^ ?° me i arge an Sle-of-attack maneuvers, merely estimating bias errors and 
scale factors for the air data may not be sufficient. In a paper d c T 

Tott^T °f in f Cial functions ’ Gupta and Iliff (ref. 37) found it necesslry 
. a . in estimates for air-variable time histories for the high angle-of-attack 
flight-test regime. The data used in the solution of this problem consisted of 
onboard measurements of Euler angles, as well as radar tracking data (slant range 
bearing, and elevation angles). Winds were estimated during low angle^-attaSk 

1 : IZITV" data “ ere USable - The ™ ~ n 

rZ IT, ► segments when the air variables were to be estimated 

estimates were obtained by "smoothing" the radar data for the Earth-frame 
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ESTIMATION OF EULER ANGLES 

ooooo TRUE SOLUTION 


ESTIMATE 



TIME, sec 


Figure 4.2.- Angles for example 2. (a) Boll angle, (b) pitch angle, 

° (c) yaw angle. 


components of aircraft velocity, subtracting the winds 
aircraft body-frame system to calculate the desired es 


, and then transforming to the 
timates of airspeed, angle of 


attack, and angle of sideslip. 


The measurement 
variables, shown in 


set employed by Gupta and II iff for 
table 4.1, is the basis for example 


the 
3. I 


estimation of air-data 
n the solution of this 
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example, the wind components are again assumed to be known. Here the estimated part 
of the forcing-function vector consisted of elements (<j>, 0, ip) and (x, y, h), and 
the measurements fitted were (<t>, 0, ip) and (R, B, h). The results of the solution 
for the air variables (V, a, 6) are shown in figure 4.3 and table 4.2. In an appli- 
cation such as this, the radar data-sample rate may not be high enough to provide 
sufficient air-variable estimates (in the Gupta application, the sample rate was 
1 Hz). It may be both necessary and convenient to augment the measurement set with 
onboard accelerometer data. Example 4 illustrates this case by including (a , a , 
a z ) in the measurement set to be fitted. The results are shown in figure 4.4 anci 
table 4.2, where a comparison can be made with the results of the preceding example. 


Example 5 

The Taylor application requires air-data measurements, whereas the Gupta-Iliff 
application requires Euler-angle measurements. It would be useful in some extreme 
flight-test situations to be able to estimate both sets of variables. That this can 
be accomplished by state estimation is illustrated by a final example (example 5). 

As indicated in table 4.1, this procedure utilizes radar position data (including 
altitude), and measurements of the "strap-down" variables (linear accelerations and 
angular velocities). Results of a simulation experiment as obtained by SMACK are 
shown in figure 4.5, where good correspondence between estimated and true time- 
histories can be observed. A comparison of the estimation accuracy obtained here 
with the results of the three previous experiments can be seen in table 4.2. 
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Figure 4. 


ooooo TRUE SOLUTION 
ESTIMATE 




Air data for example 3. (a) True airspeed, (b) angle of attack, 

(c) sideslip angle. 
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25 


YAW ANGLE, deg PITCH ANGLE, deg ROLL ANGLE, deg 


OOOOO TRUE SOLUTION 
ESTIMATE 

20 

10 


o.* — r - 4 

\ ! ! / : 



TIME, sec 


OOOOO TRUE SOLUTION 
ESTIMATE 


190 



7.50 


I* 6.25 

yf 

u 

< 5.00 

I- 

< 

0 3.75 

Lii 

-I 

c D 

Z 2 50 
< 


1.25 


3 


oj o 

0) ^ 




100.0 




12.5 25.0 37.5 50.0 62.5 

TIME, sec 


75.0 87.5 100.0 


Figure 4.5.- Euler angles and air data variables for example 5. (a) Roll angle, 
(b) pitch angle, (c) yaw angle, (d) true airspeed, (e) angle of attack, 

(f) sideslip angle. 
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5. CODING A PROBLEM FOR SMACK 


The algorithm and rigid-body model utilized by SMACK, and several state estima- 
tion applications have been discussed in previous chapters. In this chapter, the 
coding procedure for analyzing flight-test and accident data will be presented. All 
coding is prepared in 80-column statement lines. The FORTRAN 77 conventions for 
integer constants (15 format), decimal constants (F10.0 format), and Hollerith 
characters (A3 format) are used. With the exception of the first, all statement 
lines in the coding list are similarly formatted in columns (col.), as follows: 

• 1 - 3 Variable name or function descriptor 

Col. 6-10 I, a r ight- justif ied integer constant 

Col. 11-15 J, a right- justif ied integer constant 

Col. 16 - 20 K, a r ight- justif ied integer constant 

Col. 21 - 25 L, a right-justified integer constant 

Col. 31 - 40 VAL1, a decimal constant 

Col. 41-50 VAL2, a decimal constant 

Col. 51 - 60 VAL3, a decimal constant 

Col. 61-70 VAL4 , a decimal constant 

Col. 71 - 80, VAL5, a decimal constant 

The first statement line of the SMACK coding list must be a problem 
description, such as 

CODING LIST FOR AN A/C STATE ESTIMATION PROBLEM 

a message that may contain 48 characters. The next line in the list must be a 
solution description, which should be coded as 

123 10 15 20 25 

MKS I J K L 

or 

ENG I J K L 

where MKS and ENG define the system of units used to display problem variables. 
Integers I, J, K, and L are interpreted as follows: 

I number of iterations to obtain a final solution 

J number of iterations of a starting solution 

K=1 output format for accident analysis 

L-1 aircraft simulation (rising, 180° turn in wind) 
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Convergence of the SMACK algorithm is usually accomplished within ten iterations 
( I< 10) . The starting set of initial conditions and forcing functions determined by 
the program directly influences convergence properties. Occasionally it may be 
helpful to iterate the starting solution, which is done by using state measurements 
(when available) to evaluate the Jacobian matrices. The number of iterations 
desired is specified by the value of J. A special output format for displaying the 
results of an accident analysis is chosen by setting K=1. A test problem included 
in appendix D illustrates this option. The last parameter (L=1) is used to initiate 
an analysis of a simulated maneuver, which is useful for testing a given coding 
list. All of the examples shown in the previous chapter were prepared by using the 
simulated maneuver. 

Other statements in the coding list may be placed in any order, except for the 
END statement, which must be last. It appears as 


123 

10 

15 

END 

I 

J 

where, for 


I = -1 


page plot of starting solution 

= 1 


page plot of final solution 

--2 

,-3 

x-y plot of starting solution 

=2, 

3 

x-y plot of final solution 


J=1 analysis of coding list printed 

Plots include all output variables mentioned in the coding list. In the page plot, 
the time variable runs lengthwise on a printer page, and may continue for several 
pages The page-plot routine exists as a SMACK subroutine; the x-y plot routines 
require IMSL and DISSPLA libraries. The IMSL plots (I=-2,2) are produced on the 
system line printer, whereas DISSPLA plots (I=-3,3) are produced on a plotting 
device. Specification of a coding list analysis (J=1) is useful for detecting 
coding errors and as an aid in learning how the program works. 

The coding list must have one entry for each quantity considered as an "output" 
variable in the solution. Outputs include measured variables as well as variables 
to be estimated. Specification of an output variable appears as 

123 10 15 20 25 31-40 41-50 51-60 61-70 71-80 

VAR I J K L VAL1 VAL2 VAL3 VAL4 VAL5 

where 
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VAR 

variable 

name, chosen from list of table 5.1 

1 = 1 

variable 

has been measured 

J=1 

estimate 

variable time history 

K= 1 

estimate 

instrument bias error 

L = 1 

estimate 

instrument scale factor 


TABLE 5.1.- LIST OF OUTPUT VARIABLES 


Row 

Symbol 

Description 

Internal 

units 3 

External 

units 

MKS 

ENG 

1 

PHI 

Roil angle 

r 

d 

d 

2 

THT 

Pitch angle 

r 

d 

d 

3 

PS I 

Yaw angle 

r 

d 

d 

4 

X 

Position (north) 

m 

nm 

nm 

5 

Y 

Position (east) 

m 

nm 

nm 

6 

H 

Altitude (ASL) 

m 

m 

f 

7 

RNG 

Slant range 

m 

nm 

nm 

8 

BRG 

Bearing angle 

r 

d 

d 

9 

ELV 

Elevation angle 

r 

d 

d 

10 

WXY 

Horizontal wind speed 

m/s 

kt 

kt 

1 1 

WHD 

Horizontal wind heading 

r 

d 

d 

12 

VWD 

Vertical wind speed 

m/s 

m/s 

r/s 

13 

VT 

True airspeed 

m/s 

kt 

kt 

14 

AV 

Angle of attack 

r 

d 

d 

15 

BV 

Sideslip angle 

r 

d 

d 

16 

AX 

Body specific force 

m/s 2 

g 

g 

g 

g 

d/s 

17 

AY 

Body specific force 

m/s 2 

O 

g 

18 

AZ 

Body specific force 

m/s 2 

o 

g 

19 

P 

Roll rate 

r/s 

d/s 

20 

Q 

Pitch rate 

r/s 

d/s 

d/s 

21 

R 

Yaw rate 

r/s 

d/s 

d/s 

22 

AL 

Roll acceleration 

r/s 2 

d/s 2 

d/s 2 

23 

AM 

Pitch acceleration 

r/s 2 

d/s 2 

d/s 2 

2H 

AN 

Yaw acceleration 

r/s 2 

d/s 2 

d/s 2 

25 

RNA 

Slant range (aux) 

m 

nm 

nm 

26 

BRA 

Bearing angle (aux) 

r 

d 

d 

27 

ELA 

Elevation angle (aux) 

r 

d 

d 

28 

HDG 

Heading angle 

r 

d 

d 

29 

VGR 

Groundspeed 

m/s 

kt 

kt 

30 

TRK 

Groundtrack 

r 

d 

d 

Flight data must 
Appendix E). 

■ be converted to internal 

units in subroutine DATA 

(see 
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VAL1 one-sigma value of measurement noise 

VAL2 a priori estimate of instrument bias error 

VAL3 one-sigma value for a priori bias estimate 

VAL4 a priori estimate of instrument scale factor 

VAL5 one-sigma value for an a priori scale-factor estimate 

The units for VAL 1 , VAL2 , and VAL3 are designated by the choice of MKS or ENG (see 
table 5 1) The default value for either VAL3 or VAL5 is infinity (indicating no 
confidence' in the a priori estimate). Note that for external data records, setting 
VAL1 >0 takes precedence over a program-determined noise RMS value. For simulated 
(internally generated) data, setting VAL1>0 specifies the amount of noise to be 
added to the record, as well as the weight to be used in the solution. For VAL1=0, 
no noise is added and a weight of unity is used in the solution. 

Note that the preceding statement-line format should also be used when it is 
desired to use measured linear accelerations and/or angular velocities as forcing 
functions. This may be useful when performing an initial data-consistency check. 
Each measured forcing function should be specified by setting 1=1, J=0 in the 
coding list. Bias errors and scale factors may also be specified. However, it is 
usually desirable to estimate all forcing-function time histories. An entry in the 
coding list for each forcing function to be estimated should appear as 


123 10 15 20 

31-40 

41-50 

51-60 

VAR 2 J K 

VAL1 

VAL2 

VAL3 


where 


VAR 

variable 

name, chosen from list of table 

5.2 

J=1 

estimate 

forcing-function time history 


K= 1 

estimate 

forcing-function mean value 


VAL1 

RMS value of forcing function 


VAL2 

a priori 

estimate of forcing-function mean value 

VAL3 

one-sigma value for a priori mean value 

estimate 


The units for VAL1 , VAL2 , and VAL3 are designated by the choice of MKS or ENG. The 
default value for VAL3 is infinity. Note that if VAL1>0, that value takes prece- 
dence over a program-determined RMS value for either external or simulated data 
records. When all the elements of either (AL, AM, AN) or (AX, AY, AZ) are specified 
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TABLE 5.2.- LIST OF FORCING-FUNCTION VARIABLES 


Row Symbol Description 


Internal External 
units units 

MKS ENG 


1 

DL 

PH2 time-derivative 

r/s 3 

d/s 3 

d/s 3 

2 

DM 

TH2 time-derivative 

r/s 3 

d/s 3 

d/s 3 

3 

DN 

PS2 time-derivative 

r/s 3 

d/s 3 

d/s 3 

4 

DX 

X2 time-der ivative 

m/s 3 

m/s 3 

f/s 3 

5 

DY 

Y2 time-derivative 

m/s 3 

m/s;: 

f/s 3 

6 

DH 

H2 time-derivative 

m/s 3 

m/s 3 

f/s 3 

7 

GX 

WX time-derivative 

m/s 2 

m/s 2 

f/s 2 

8 

GY 

WY time-derivative 

m/s 2 

m/s 2 

f/s 2 

9 

GH 

WH time-derivative 

m/s 2 

m/s 2 

f/s 2 

10 

PH2 

PHI time-derivative 

r/s 2 

d/s 2 

d/s 2 

11 

TH2 

TH1 time-derivative 

r/s 2 

d/s 2 

p 

d/s 

12 

PS2 

PS1 time-der ivative 

r/ s 2 

d/s 2 

d/s 2 

13 

X2 

X 1 time-der ivative 

m/s 2 

m/s 2 

f/s 2 

14 

Y2 

Y 1 time-derivative 

m/s 2 

m/s 2 

f/s 2 

15 

H2 

H 1 time-der ivative 

m/s 2 

m/s 2 

f/s 2 


with 1=1, J=1 in the coding list, forcing functions (DL, DM, DN) or (DX, DY, DH) are 
selected by the program. Hence, those forcing functions should be included in the 
list only to override a program-determined RMS weight, or to specify estimation of a 
mean value. Their inclusion in the coding list, however, is a useful reminder of 
the excitations chosen for the state model shown in figure 3.1. Care should be 
taken not to mix elements of (DL, DM, DN) with elements of (PH2, TH2, PS2) or 
elements of (DX, DY, DH) with elements of (X2, Y2, H2). 

Inclusion of a particular state variable in the coding list is necessary only 

to specify an a priori initial condition of that state. The statement line should 
appear as 


123 10 41-50 51-60 

VAR 3 VAL2 VAL3 

where 


VAR 

VAL2 

VAL3 


variable name, chosen from list of table 5.3 
an a priori estimate of the initial condition 
one-sigma value for the a priori state estimate 
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TABLE 5.3.- LIST OF STATE VARIABLES 


Row Symbol Description Internal External 

units units 

MKS ENG 


1 

PH2 

PHI time-derivative 

r / sz 

d/s^ 

d/s 

2 

TH2 

TH1 time-derivative 

r/s^ 

d/s^ 

d/s 

3 

PS2 

PSI time-derivative 

r/s^ 

d/s^ 

d/s 

4 

X2 

XI time-derivative 

m/ s'z 

m/s^ 

f /s 

5 

Y2 

Y1 time-derivative 

m/s^ 

m/ s'z 

f/s 

6 

H2 

HI time-derivative 

m/s^ 

m/s d 

f/s 

7 

WX 

Wind speed (north) 

m/s 

m/s 

f/s 

8 

WY 

Wind speed (east) 

m/s 

m/s 

f/s 

9 

WH 

Wind speed (up) 

m/s 

m/s 

f/s 

10 

PHI 

PHI time derivative 

r/s 

d/s 

d/s 

1 1 

TH1 

THT time-derivative 

r/s 

d/s 

d/s 

12 

PS1 

PSI time-derivative 

r/s 

d/s 

d/s 

13 

XI 

X time-derivative 

m/s 

m/s 

f/s 

14 

Y 1 

Y time-derivative 

m/s 

m/s 

f/s 

15 

HI 

H time-derivative 

m/s 

m/s 

f/s 

16 

PHI 

Roll angle 

r 

d 

d 

17 

THT 

Pitch angle 

r 

d 

d 

18 

PSI 

Yaw angle 

r 

d 

d 

19 

X 

Position (north) 

m 

nm 

nm 

20 

Y 

Position (east) 

m 

nm 

nm 

21 

H 

Altitude (ASL) 

m 

m 

f 


It should be emphasized that specifying an a priori estimate with its corresponding 
one-sigma value for any variable will bias the solution towards that estimate. 

A priori estimates will seldom be necessary for convergence of the algorithm. 

Again, the units for VAL2 and VAL3 are designated by the choice of MKS or ENG. 

Each coding list must contain a description of the data record and the way in 
which it is to be processed. This description includes the number of data points, 
the sampling interval, and the integration time step. The code should appear as 

123 10 15 20 25 31-40 41-50 

REC I J K L VAL1 VAL2 

where 

I starting point of record 

J ending point of record 
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K number of integration steps per sampling interval 

L integration steps per output point (plotting only) 

VAL1 data sampling interval in seconds 
VAL2 filter cutoff frequency in Hertz 

For measurement sets with multiple-rate data, the starting and ending points should 
correspond to the record with the highest data rate. The time step used for inte- 
gration is chosen to be an integral submultiple of each of the sampling intervals 
(see appendix E). The sampling interval VAL1 should be that of the record with the 
highest data rate. Mote that each measurement record is low-pass filtered in the 
starting routine in order to obtain a measure of the residual covariance. The 
cutoff frequency is VAL2, which should be adjusted so that residuals of the records 
with the highest data rate are as "white" as possible. The cutoff frequencies for 
other records in the measurement set are program-scaled by sample-rate ratios. The 
default value for VAL2 is 0.1/VAL1 Hz. 

It may sometimes be useful to independently specify xy-plot scales or a filter 
cutoff frequency for a record. This can be done by including in the coding list the 
statement line 

123 10 31-40 41-50 51-60 61-70 71-80 

VAR 4 VAL1 VAL2 VAL3 VAL4 VAL5 

where 

VAR variable name, chosen from list of table 5.1 

VAL1 filter cutoff frequency in Hertz 

VAL2 x-axis minimum value 

VAL3 x-axis maximum value 

VAL4 y-axis minimum value 

VAL5 y-axis maximum value 


The units for VAL2 through VAL5 are designated by the choice of MKS or EMG. If all 
four values are zero, no changes in program-determined plot scales will be made. If 
VAL1 is zero, the filter cutoff frequency for the record will be that specified by 
the REC statement line. Note that the filter residuals may be examined for white- 
ness by obtaining plots of the starting solution ( I =— 1 in the END statement). 
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The SMACK program can compensate for instrument offset from the aircraft c.g. 
One correction can be made for each of three instrument sets: accelerometer, pres- 

sure ports, and alpha-beta vanes. These statement lines should appear as 

123 41-50 51-60 61-70 

ACC VAL2 VAL3 VAL4 

for the accelerometer package, 

P-S VAL2 VAL3 VAL4 

for the pitot-static system, and 

VNE VAL2 VAL3 VAL4 

for the vane measurement system. In each statement, VAL2, VAL3, and VAL4 represent 
the (x, y, z) location of the instrument set with respect to the c.g., in meters 
(MKS) or in feet (ENG) . 

For a problem that includes radar tracking data, it is possible to specify the site 
location with respect to a desired origin using the statement line 


123 

41-50 

51-60 

61-70 

RAD 

VAL2 

VAL3 

VAL4 


where VAL2, VAL3, and VAL4 represent the (x, y, h) location of the tracking 
antenna. Here VAL2 and VAL3 are in nautical miles, and VAL4 is in meters (MKS) or 
feet (ENG). When an auxiliary site has provided tracking data, its location can 
similarly be represented with 

RDA VAL2 VAL3 VAL4 

For accident analysis, when the only data available are the radar track (including 
altitude), winds, air temperature and aircraft performance data, none of the air- 
craft trajectory variables (e.g., attitude or velocity) need be specified in the 
coding list. The trajectory variables will be determined following the radar solu- 
tion by specifying K=1 in the solution description. The first problem included in 
appendix D illustrates this application. 

Preparation of the coding list will here be illustrated by returning to the 
application examples of the previous chapter. A list for each example is found in 
figures 5.1 through 5.5. In each case, the solution description specifies MKS 
units, eight iterations and the SMACK aircraft simulation to provide the data 
records. The REC statement specifies that each record will have 90 points, the 
integration step will be the same as the sampling interval, every point will be 
plotted, and the sampling interval is one second. To help interpret the VAR 
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EXAMPLE Is ESTIMATION OF WINDS 


MKS 

8 


1 


REC 

1 

90 

1 1 

1.0 

AX 

1 

1 

1 

0.001 

AY 

1 

1 

1 

0.001 

AZ 

1 

1 

1 

0.001 

RNG 

1 

1 


0.001 

BRG 

1 

1 


0.05 

H 

1 

1 


0.5 

WXY 


1 



WHD 


1 



VWD 


1 



GX 

2 

1 



GY 

2 

1 



GH 

2 

1 



PHI 

1 

1 


0.05 

THT 

1 

1 


0.05 

PSI 

1 

1 


0.05 

PH2 

2 

1 



TH2 

2 

1 



PS2 

2 

1 



VT 

1 

1 


0.1 

AV 

1 

1 


0.05 

BV 

1 

1 


0.05 

END 

2 

1 



F igure 5.1 

.- Coding 

list 

for application 

example 


EXAMPLE 2: ESTIMATION OF EULER ANGLES 


MXS 

REC 

AX 

AY 

AZ 

P 

<3 

R 

WXY 

WHD 

VWD 

PHI 

PHI 

THT 

THT 

PSI 

PSI 

PH2 

TH2 

PS2 

VT 

AV 

BV 

END 


8 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

3 

3 

3 

2 

2 

2 

1 

1 

1 

2 


90 1 

1 
1 
1 

1 1 

1 1 

1 1 


1 

1 

1 

1 

1 

1 

1 1 
1 1 
1 1 
1 


1 

1 


1 

1 


1.0 

0.001 

0.001 

0.001 

0.05 

0.05 

0.05 


0. 

1.31 

90. 


0.1 

0.05 

0.05 


5 

5 

1 


Figure 5.2.- Coding list for application example 2. 
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EXAMPLE 3: ESTIMATION OF AIR VARIABLES 


MKS 

8 


1 


REC 

1 

90 

1 1 

1.0 

RNG 

1 

1 


0.001 

BRG 

1 

1 


0.05 

H 

1 

1 


0.5 

X2 

2 

1 

1 


Y2 

2 

1 

1 


H2 

2 

1 

1 


WXY 

1 



0.1 

WHD 

1 



0.05 

VWD 

1 



0.1 

PHI 

1 

1 


0.05 

THT 

1 

1 


0.05 

PSI 

1 

1 


0.05 

PH2 

2 

1 



TH2 

2 

1 



PS2 

2 

1 



VT 


1 



AV 


1 



BV 


1 



END 

2 

1 




Figure 5.3.- Coding list for application example 3 


EXAMPLE 4: ESTIMATION OF AIR VARIABLES 


MKS 

8 

1 


REC 

1 

90 1 1 

1.0 

AX 

1 

1 

0.001 

AY 

1 

1 

0.001 

AZ 

1 

1 

0.001 

RNG 

1 

1 

0.001 

BRG 

1 

1 

0.05 

H 

1 

1 

0.5 

WXY 

1 


0.1 

WHD 

1 


0.05 

VWD 

1 


0.1 

PHI 

1 

1 

0.05 

THT 

1 

1 

0.05 

PSI 

1 

1 

0.05 

PH2 

2 

1 


TH2 

2 

1 


PS2 

2 

1 


VT 


1 


AV 


1 


BV 


1 


END 

2 

1 


Figure 5.4. 

- Coding 

list for application exampL 
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EXAMPLE 5: ESTIMATION OF ANGLES, AIR VARIABLES 


MKS 

REC 

AX 

AY 

AZ 

P 

9 

R 

RNG 

BRG 

H 

WXY 

WHD 

VWD 

PHI 

PHI 

THT 

THT 

PSI 

PSI 

PH2 

TH2 

PS2 

VT 

BV 

AV 

END 


8 

1 90 

1 1 

1 1 

1 1 

1 1 

1 1 

1 1 

1 1 

1 1 

1 1 

1 
1 
1 



1 

3 

2 1 

2 1 

2 1 

1 
1 
1 

2 1 


1 


1 

1 

1 


1 

1 


1.0 

0.001 

0.001 

0.001 

0.05 

0.05 

0.05 

0.001 

0.05 

0.5 


0. 

2.0 

1.313 

to 

o 

91.35 

1.0 


Figure 5.5. Coding list for application example 5. 

statement lines, refer again to table 4.1, which lists the measured and/or estimated 

SeluJer flf '“"P 1 *; ***“ that in examples 2 and 5, a priori estimates of 

g es W6re specifled - For example 3, forcing functions (X2, Y2, H2) were 
specified, but in example 4, accelerometer measurements (AX, AY, AZ) were available 
and corresponding forcing functions (DX, DY, DH) were used by the program In 
example 4, the user does not need to specify the forcing-function set For 

es^L 1 ^ 2 T°!f 5 ’ ^ WindS (WXY ’ WHD ’ VWD) Were ^ified as measured, but not 
estimated. In these examples, states (WX, WY, WH) are computed from the measure- 

ments and used with (XI, Y1, HI) in the estimation of air velocities. Finally note 

a the coding list for each example concludes with an END statement specifying 

Te to'Z printed!" °“ tPUt VariableS and a diagnostic analysis of the problem 
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6. A FLIGHT-TEST METHODOLOGY 


This chapter describes a flight-test methodology for acquiring a data base to 
identify a full-envelope aerodynamic model of a V/STOL Research Aircraft (VSRA). 

The model will serve to update and improve an existing VSRA simulation, in order to 
aid the design of guidance, control, and display systems for the aircraft. A key 
element in the methodology is the application of SMACK for the processing of each 
test maneuver before its entry to the data base. It should be helpful to the reader 
to see how the state-estimation method may be used in a flight-test setting. For a 
more complete discussion see references 38 and 39- 

The NASA VSRA is a YAV-8B aircraft, a prototype of the subsonic, vectored- 
thrust "Harrier" fighter aircraft; its engine nozzles can be rotated from zero 
degrees for forward flight to somewhat greater than 90° for hover and vertical 
flight. A reaction-control system (RCS), in which compressor air is piped to the 
extremities of the aircraft, provides attitude control in hover and low-speed 
flight. The VSRA aerodynamic model must represent the three body forces and three 
moments over a flight envelope that includes hover, transition to forward flight and 
back to hover, and STOL operation and normal cruise. 

The resulting model, strongly nonlinear with respect to aircraft variables such 
as angles of attack and sideslip, Mach number, nozzle angle, and power setting, can 
be conveniently expressed with functions that are linear in the parameters to be 
identified (refs. 40-42). A linear least-squares (regression) method (refs. 43-4b) 
is well-suited to identify a highly nonlinear model that is linearly parame- 
terized. Because regression methods are computationally simple, careful attention 
can be given to the structuring of an accurate and physically meaningful model. 

Good results with regression methods, however, are highly dependent on the quality 
of the flight data. Therefore, state-estimation methods are often used before 
modelling to correct the data records for bias and scale-factor errors and to pro- 
vide estimates of unmeasured or poorly measured variables. 


The methodology for acquiring a data base matched to a least-squares (regres- 
sion) identification task is outlined in the flow diagram shown in figure 6.1. The 
important aspects of the preflight planning, flight testing, and postflight Press- 
ing phases necessary to acquire the data base are covered in this chapter. e ro e 
of state estimation in the processing is emphasized. The actual modelling of V 
aerodynamics is beyond the scope of this discussion. 


Maneuver Design 

The data base required for least-squares aerodynamic model identification can 
be obtained quite efficiently. Because the model is nonlinear, it is not necessary 
(or useful) to maintain trim during a maneuver. In addition, because a regression 
procedure will be used to identify the model, large amounts of data may be batch- 
processed. Accordingly, each flight-test maneuver has been designed to yield large 
changes in aircraft variables while covering a (nearly) closed course within five 
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RADAR 

TRACKING 




Figure 6.1.- Flow diagram for VSRA flight-test methodology. 


39 











minutes under continuous radar tracking. The raw data base consists of as many 
longitudinal, lateral, and transitional (to and from hover) maneuvers as are neces- 
sary to cover the flight envelope. After processing, model sections may be iden- 
tified using long (15-30 min) records, each consisting of concatenated segments from 
several maneuvers. 

One set of maneuvers was designed to excite large changes in longitudinal-model 
variables from several nominal trim points. In each of these maneuvers, the nozzle 
angle, flap deflection, and power are held constant while the stabilator is varied 
to obtain changes in angle of attack and pitch rate. The maneuver includes "stick 
pulses," sinusoidal "stick pumping," an "alpha ladder," and a "wind-up" turn 
(ref. 39). Wear the end of the maneuver, power is added to return to the nominal 
trim point. Note that a considerable variation in Mach number may be experienced 
during the maneuver. Another set of maneuvers was designed to excite large changes 
in lateral-model variables (angle of sideslip, yaw rate, and roll rate). Most of 
the maneuvers were performed without "stability augmentation" to ensure a full range 
of aircraft response activity. All V/STOL procedures were performed in and out of 
ground effect. 

One characteristic that sets the VSRA apart from conventional aircraft is that 
it exhibits significant thrust-induced aerodynamic effects when the nozzles are not 
in the full-aft position. These are largest during transition from hover to forward 
flight (and back to hover) and during periods of low-speed flight. Standard V/STOL 
procedures were used to provide data for identification of thrust-induced aerody- 
namics. One of these procedures, a short-takeoff and slow-landing maneuver, is out- 
lined on the flight-test card shown in figure 6.2. In this maneuver, the ground 
roll begins with nozzles at 10°. At V p (indicated air speed) the nozzles are 
rotated to an angle 0 p (in the example for this chapter, V r = 50 kt and 
9 r = 55°). Shortly after liftoff, the nozzles are rotated to the full-aft posi- 
tion. For the slow-landing portion, nozzles are rotated to 40° just before the 
final turn, and during the final approach are further rotated to 60°. 


Data Acquisition 

The VSRA measurement system is equipped with a 10-bit digital data acquisition 
and telemetry (TM) system. A pulse-code modulation format is used to encode 156 
mainframe channels sampled at 120 Hz and 160 subframe channels sampled at 30 Hz. 
Before encoding, each analog channel is passed through a third-order Butterworth 
anti-aliasing filter with its cutoff frequency set at one-fifth of the channel 
sampling rate. After encoding, all flight data are transmitted to a ground station 
where they are recorded. A partial list of onboard measurements, those necessary 
for aerodynamic model identification, is given in table 6.1. 

Flight tests of the VSRA were performed at the NASA test facility located at 
Crows Landing, California. The facility control room, which has a clear view of the 
runway and hover pad, is equipped with five eight-channel strip-chart recorders and 
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FLIGHT TEST CARD 

SHORT-TAKEOFF SLOW-LANDING 

Aircraft: VSRA (NASA 704) Flight: 744 

Experimenters: McNally/Bach Date: 11/12/87 


Pilot: Gerdes 



Figure 6.2.- Plan view of example VSRA flight-test maneuver (exact reproduction 

of flight-test card used by pilot). 
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TABLE 6.1.- VARIABLE LIST FOR AERODYNAMIC MODEL DATA BASE 


Channel 

Measured 

Estimated 

Euler angles 

Onboard 

SMACK 

Angular rates 

Onboard 

SMACK 

Angular accelerations 


SMACK 

Linear accelerations 

Onboard 

SMACK 

Inertial positions 

Radar 

SMACK 

Inertial velocities 


SMACK 

Air-flow angles 

Onboard 

SMACK 

Static pressure 

Onboard 


Total pressure 

Onboard 


Total temperature 

Onboard 


True airspeed 


SMACK d 

Flightpath winds 


SMACK 

Flap setting 

Onboard 


Aileron deflections 

Onboard 


Stabilator deflection 

Onboard 


Rudder deflection 

Onboard 


Engine nozzle angle 

Onboard 


Engine fan speed 

Onboard 


Compressor pressure 

Onboard 


Fuel and water weights 

Onboard 


RCS roll-valve positions 

Onboard 


RCS pitch-valve positions 

Onboard 


RCS yaw-valve position 

Onboard 


Engine and RCS body forces 


ENCAL 

Engine and RCS moments 


ENCAL 

Gross weight and inertias 


ENCAL 


a SMACK utilizes a "measurement" of true airspeed, which 
is computed from the ratio of totax and static pressures, 
and the total temperature (see chapter 3 and ref. 29). 

three color monitors for real-time display of the TM data. Two on-site radar 
systems are available to provide continuous tracking of the test aircraft posi- 
tion. (A laser tracking system is used for all hover maneuvers.) During flight 
test, TM data from the VSRA onboard system are downlinked, merged at the facility 
with range, bearing, and elevation data from the tracking systems, and then 

recorded . 
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Preliminary Processing 


Following real-time acquisition of data during flight test, the records from 
each maneuver are converted to engineering units and made available to researchers 
in a raw flight-data file. The first program in the postflight processing system 
reads the raw file and creates a "processed" file of selected channels. The pro- 
cessing begins by removing wild points from the records. Several options are avail- 
able, but one effective (but time-consuming) method is to pass each record through a 
"moving window". Points that fall outside the window are considered wild, and are 
tagged but not removed. When all wild points in a record have been tagged, the 
record is passed through a low-pass digital filter (see appendix E) to obtain an 
interpolated time history free of wild points. After interpolation, the data rate 
can be reduced to a submultiple of the mainframe sampling frequency. The filter 
cutoff frequency is set at one-half the final data rate desired. The final rate was 
chosen to be 20 Hz for all VSRA maneuvers. 

Each channel processed from a maneuver raw-data file is stored in a processed 
flight-data file set up for that maneuver. The analyst now may use a program to 
interactively select processed data channels for plotting in either x-y or strip- 
chart format. An x-y cross plot, for example, might display Mach number plotted 
against angle of attack. Such plots offer a convenient way to evaluate how well the 
flight envelope has been covered during a maneuver. It is unlikely that a single 
maneuver will provide enough variation in aircraft variables to identify all model 
terms: the analyst may also use this program to create a "map" file, which will 

contain addresses of time segments selected from several processed maneuvers. This 
file can later be used to concatenate the selected segments to create a long record 
suitable for model identification. 


State Estimation 

The next step in the processing of each maneuver is to apply SMACK to check 
data consistency and derive unmeasured variables from the measurement set given in 
table 6.1. The relatively long (3-5 min) maneuvers with large dynamic variations 
are well-suited to state-estimation analysis. The closed course yields good track- 
ing accuracy and facilitates estimation of winds along the flightpath. The aircraft 
dynamic response to the control inputs is analyzed (along with the radar track) by 
the SMACK procedure, which determines integrator initial conditions, selected 
instrument bias errors and scale factors, and forcing-function time histories that 
provide the "best fits" to the measurement records. The body angular accelerations, 
true airspeed, and flightpath winds are also estimated as part of the solution. 

The user must write and link subroutine DATA with the SMACK program. This 
subroutine accesses the data records and performs the chores required to prepare the 
problem for solution. In addition to setting up the arrays according to the list in 
table 5.1, for example, there may be intervals during which some instruments are 
known to be "saturated" or else sections of data may be unusable for some other 
reason. The data in such intervals must be "blanked" so that they will not 
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influence the estimates. These chores are relatively easy to accomplish; the 
procedures are described in appendix E. 

A coding list for the analysis of a typical VSRA flight-test maneuver is shown 
in figure 6.3. The REC statement indicates a data record of 6000 points and a 
sampling interval of 0.05 sec (these may be changed by the user in subroutine 
DATA). A filter cutoff frequency of 6 Hz has been specified for determination of 
performance- index weights. Notice that data from both tracking systems are being 
used; the RAD and RDA statements specify their positions with respect to a runway 
origin. Notice also that position corrections for accelerometer and air-data 
instrument locations on the aircraft are specified by the ACC, P-S, and VNE 
statements . 


Force and Moment Calculations 

The aerodynamic forces and moments acting on the VSRA during flight are deter- 
mined as the difference between the total forces and moments and the engine forces 
and moments. Here the term "engine" includes the reaction control system as well as 
the main nozzles. The engine forces and moments are calculated (offline) by a 
program called ENCAL (ENgine CALculations) . This program uses a nominal propulsion 
model of the VSRA Pegasus engine (YF402-RR-404) (ref. 46). Fan dynamics are not 
included in this version, since fan speed is measured in flight. It should be noted 
that the propulsion model provides only thrust forces and moments. Any thrust- 
induced aerodynamic effects are to be included in the VSRA aerodynamic model. 

Inputs to the ENCAL routine include all the air-data, reaction-control, engine, 
and weight measurements listed in table 6.1. Outputs to the processed flight-data 
file are the three body-axis components of engine force and moment. The ENCAL 
routine also calculates aircraft weight and inertias, and the variation in center- 
of-gravity location. These variables are added to the processed-data file. Note 
that the aerodynamic model to be identified from flight data can only be as accurate 
as the engine model. A fully-instrumented Pegasus engine has recently been 
installed on the VSRA, and the engine model will be validated after the next set of 
flight tests. 

Total VSRA force and moment time-histories are obtained from the SMACK-derived 
estimates of accelerations and angular rates, and from ENCAL-der ived estimates of 
weight and inertias. The body-axis forces are given by 

F x = ma x ; F y = ma y ; F z = ma z (6.1) 

where m is vehicle mass. The moments are calculated from 


T = 

*xx a Sl 

^zx^ a n 
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(a) VSRA FLIGHT-TEST ANALYSIS 


(INPUTS ESTIMATED) 


ENG 

REC 

PHI 

THT 

PSI 

P 

Q 

R 

AL 

AM 

AN 

DL 

DM 

DN 

RNG 

BRG 

ELV 

RAD 

RNA 

BRA 

ELA 

RDA 

AX 

AY 

AZ 

ACC 

DX 

DY 

DH 

END 


10 0 0 0 
1 6000 1 1 

110 1 
110 0 

110 0 

1110 
1110 
1110 
0 10 0 

0 10 0 

0 10 0 

2 10 0 

2 10 0 

2 10 0 

110 0 

110 0 

110 0 

oooo 

110 0 

110 0 

110 0 

oooo 
1110 
1110 
1110 
oooo 
2 10 0 

2 10 0 

2 10 0 

3 10 0 


0.05 

6. 

0.2 

0.8 

0.2 

-0.6 

0.5 

-2.3 

0.1 


0.1 


0.1 



0.005 

0.05 

0.05 

0 • .28337019 .19029757 

0.005 

0.05 

0.05 

°* .33024389 .14118231 

0.002 

0.002 

0.005 

0. 0.58 0. 


(b) 


VT 

P-S 

AV 

BV 

VNE 

WXY 

WHD 

VWD 

GX 

GY 

GH 


110 0 

0000 
1110 

1 1 1 o 

oooo 

0 10 0 

0 10 0 

0 10 0 

2 10 0 

2 10 0 

2 10 0 


0.1 

0. 25. 0. 

0.1 0 . 0 . 

0.1 0 . 0 . 

0. 25. o. 


Figure 6.3.- Coding list for VSRA data-consistency analysis, (a) Inertial 
only, (b) additional statements for full solution. 


1 . 59 


1 . 58 


2.58 


1.670 
1 . 136 
1 . 333 
1.670 


solution 
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where 


ttti are vehicle moments of inertia, 
xx’ yy’ zz* zx 


Example Maneuver 

The short-takeoff and slow-landing maneuver described earlier illustrates the 
type of information that is stored in the VSRA data base. The maneuver contains 
abrupt changes in nozzle and flap angles. The aircraft transitions to normal flight 
after takeoff, performs a "go-around," and then transitions back to a STOL configu- 
ration for a slow landing. The raw data file includes all the onboard inertial and 
air data, and radar tracking measurements as indicated in table 6.1. The variations 
in nozzle angle, flap setting, power, and control-surface positions required to 
perform the maneuver are shown in figure 6.4. Only the left aileron is shown: both 

ailerons are set to 15° down (drooped) during takeoff and landing. Note how these 
time histories correlate with the activity requested of the pilot on the flight-test 
card of figure 6.2. 

Results of the SMACK analysis required for calculating forces and moments are 
shown in figure 6.5. During a preliminary solution, a large error was noticed in 
the fit of longitudinal acceleration (AX) during the takeoff portion of the maneu- 
ver. The accelerometer had saturated at 0.6 g's, and its output in that interval 
had to be "blanked". Fortunately, the good tracking data provided the redundancy 
necessary to yield the estimate during the blanked interval. The other fits to the 
measurement time histories were quite good. Although there are no measurements of 
angular accelerations in the processed data file, the measurements of angular rates 
are of sufficient quality to ensure confidence in the acceleration estimates. 

It should be noted that the large activity in the angular accelerations of 
figure 6.5(b) is related to the reduced damping of the aircraft without stability 
augmentation. In effect, the pilot must provide the control inputs to stabilize the 
aircraft. The control-surface motions in these test data are well-correlated with 
the angular accelerations. A similar maneuver flown with stability augmentation 
shows significantly smaller excursions. The larger aircraft response activity 
obtained without augmentation will, of course, enhance the "identifiability of the 
aerodynamic model. 

As a final step in maneuver processing, the aerodynamic forces and moments are 
calculated as the difference of total and engine forces and moments as outlined in 
the previous section. These are the time histories that must be adequately repre- 
sented by the VSRA aerodynamic model. Results of the ENCAL calculations for the 
maneuver are shown in figure 6.6, with the corresponding aerodynamic variables shown 
in figure 6.7. Notice the tradeoff between engine and aerodynamic vertical forces 
during the STOL portions of the maneuver. 
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Figure 6.4.- Control variations, VSRA flight-test maneuver, (a) Nozzle, flap, 

power, (b) stabilator, aileron, rudder. 
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7 . THE ANALYSIS OF WINDSHEAR 


™ ls chapter expands on the application of state estimation for the analysis of 
windshear, which was the subject of the first example of chapter 4 Encounters with 

undent H IT* represent a c ™“"hing safety problem that must be better 

bj antw™ mghfd«a diStU t b tT S that aff6Ct alrUne operatI ° ns be studied 
y g g fca recorded during typical encounters. In the past such 

usuaUy tampered by the lack of good data, but more recent 

recorder,^ lnvolved equipped with digital flight-data 

( atp ) raH ", DR records » together with ground-based air-traffic control 

from flight UsT tTST"?*' aPPr ° aChi "® that 

lence, and to characterise ^ , tS^^S.“ rCr * ft P,rf "“ , “ la 

In assisting the National Transportation Safety Board fWT^R’i in it- • <. • 

tions of accidents involving aircraft not equipped 

path “ Center developed methods to determine aircraft motions along a flight- 
path from the limited data available following an accident (refs. 33 34 and 47) 

r=rr°"lnrr b % d f e ™ ined “ ltl > SMACK state-estiLtiin method' 

. _ . ^ ). In studies of turbulence encounters involving DFDR-eauinnpri 

airiners, SMACK has been applied to determine winds along the ?iightpa?h A oar 

given In taLeTr 8 ‘ n ° lde " tS that A " es has “alymed with theOTSB is 

ence ,CA t ? i ^t^^^ 

windshear VVTX °\ 

typical CAT encounter is illustrated in figure 7 . 7 ( 1 ) ’ ' 9K # 

A more hazardous type of atmospheric disturbance is a "downburst " which i« a 
strong, concentrated downflow that induces a high-velocity outflow wirh • hll? 

ZtTr'M Z eTT- fl d typioal imbedded 

Airlines Flight £ ^ 

case an L-i 0" on final approach flew Into a downburst and 8 "^ evicted the 

encounter 3 "'!^ ^is S organized U as r foliows^ Sm ^e d n^ d i 8y ”| Pbhp,laapdi ^ht d8 9 8r downburst 

merging and synchronizing the flight data; su^quenf ^t?“ 3 r ^llnTtte°n”-dr f a 

in r t b hi e i n a°s e ; s^ir ltd ° rtb '" F ““ 
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TABLE 


Figure 7.1. 


7 . 1 .- AIRLINE TURBULENCE ENCOUNTERS REPORTED TO THE NTSB 
AND INVESTIGATED AT AMES RESEARCH CENTER 


Case 

Aircraft 

Location 

Date 

1 

DC-10 

Hannibal, MO 

4/81 

2 

DC-10 

Morton, WY 

7/82 

3 

DC- 10 

Near Bermuda 

10/83 

4 

L-101 1 

Offshore SC 

11/83 

5 

DC- 10 

Calgary, AL 

11/75 

6 

B-747 

Over Greenland 

1/85 

7 

B-747 

Over Greenland 

2/85 

8 

B-747SP 

Offshore CA 

2/85 

9 

L-101 1 

Dallas/Ft. Worth 

8/85 




- Windshear disturbances, (a) High-altitude vortex, (b) low level 

downburst . 
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Flight Data Processing 

The procedure used to determine winds along the flightpath is the same whether 
the turbulence encountered is associated with vortices at cruise altitude or with a 
downburst near the ground. In this procedure, which is outlined in figure 7.2 data 
from the DFDR and ATC radar records are synchronized and merged and the air data are 
corrected. Performance calculations necessary to synthesize unmeasured time his- 
tories of angles of attack and sideslip are made. When a complete data set is in 
place, the SMACK state-estimation algorithm is employed to estimate the winds. This 
section describes the processing required to prepare the data set for performance 
calculations and wind estimation. 


Data from a DFDR usually include measurements of accelerations, Euler angles 
pressure altitude, airspeed, and other variables, typically sampled at intervals of 
.2t>-4.0 sec. The important recorder parameters for an L-1011 DFDR data system are 
shown in table 7.2. The "frame- duration for the system is four seconds; there are 



Figure 7.2.- Estimation of flightpath winds from flight records. 
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TABLE 7.2.- PARAMETERS IN THE L- 101 1 DFDR SYSTEM 


Record 

Rate, 

Hz 

Skew 

Subframe 

Vertical acceleration 

4 

12 

1-4 

Lateral acceleration 

4 

14 

1-4 

Longitudinal acceleration 

4 

1 

1-4 

Roll angle 

1 

16 

1-4 

Pitch angle 

1 

50 

1-4 

Heading angle 

1 

2 

1-4 

Indicated airspeed 

1 

18 

1-4 

Angle of attack (1 vane) 

2 

10 

1-4 

Angle of attack (r vane) 

2 

24 

1-4 

Pressure altitude 

1 

4 

1-4 

Air temperature 

1/2 

54 

2,4 

Stabilator deflection 

1 

39 

1-4 

Rudder deflection 

2 

26 

1-4 

Thrust (engine 1) 

1/4 

32 

1 

Thrust (engine 2) 

1/4 

32 

2 

Thrust (engine 3) 

1/4 

32 

3 


four "subframes" in each frame, and each subframe has 64 sampling "slots. The 
column headed by "Rate" defines the basic sampling rate for each parameter; the one 
headed by "Skew" defines the delay (in 64ths of a second) from the start of a sub- 
frame until the parameter is first sampled. The last column specifies the sub- 
frame(s) in which the sample appears. For example, a parameter sampled at a rate of 
4 Hz with a skew of 14 would occupy slots 14, 30, 46, and 62 in each subframe. How- 
ever, a parameter sampled at a rate of 0.25 Hz with a skew of 32 would occupy slot 
32 in only one subframe of each frame. 

The first step in processing the DFDR data is to interpolate each measured 
parameter at the highest sampling rate (usually 4 Hz) before performing air-data 
corrections and other calculations. The interpolation is accomplished with a digi- 
tal filtering algorithm (see appendix E) operating at a rate of 64 Hz, in order to 
properly accommodate parameter skews. The filter also provides Euler-angle time- 
derivative estimates for use in computing body angular rates needed for estimating 
angles of attack and sideslip (or for correcting vane angles). After filtering, 
each parameter is down-sampled from 64 Hz to the appropriate rate (4 Hz) and the 
aforementioned calculations are performed. 

The second step is to correctly merge the DFDR data with the ATC radar data. 
Although each data source is time-tagged, there may be an absolute timing error of 
several seconds on either (or both) of the sources. However, there is usually 
included with the radar track an independent (transponded) record of the aircraft 
pressure altitude which can be compared with the DFDR altitude record for time- 
synchronization of the sources. Since the encoding altimeter for the transponder 
registers in increments of 100 ft, a fairly large change in altitude is necessary 
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for synchronizing the sparsely-sampled radar track with the one-Hz DFDR altimeter 
measurement record. 

The last step in data processing prior to performance calculations or wind 
estimation is to make the usual air-data computations (ref. 29). These include 
calculation of Mach number, dynamic pressure, true airspeed, and correction of the 
vane angle measurement for upwash and pitch rate to obtain the angle of attack (when 
the vane angle is included with the DFDR records). It should be noted that the 
angle-of-attack time history is essential in determining vertical wind in an inves- 
tigation of a severe turbulence encounter. 


Performance Calculations 

The time histories of force coefficients derived from flight data can be quite 
useful in accident investigations. The lift coefficient can be employed to estimate 
angle of attack (a) when that record is not among DFDR measurements (ref. 53); a 
similar procedure is generally used to reconstruct the sideslip angle (b) from a 
time history of the side-force coefficient. Both lift and drag are used in studies 
of possible performance degradation which might be caused by heavy rain or ice. 

This section reviews the use of performance calculations in analyzing turbulence 
encounters . 

Aircraft force coefficients can be expressed in two ways. In the first set, 

the lift, drag, and side-force coefficients are given in terms of measurements of 

body-axis accelerations (a x , a , a ) and thrust components (T Y , T, T ) by 

* x y z 


C L = 

[(ma x - T x )sina - (ma z - T z )cosa]/QS 

(7.1) 

C D = 

-[F cose + (ma y - T y )sinB]/QS 

(7.2) 

o 

a 

ii 

[F sine - (ma y - T y )cosB]/QS 

(7.3) 


where 


F = (ma^ - T x )cosct + (ma z - T z )sino 

where m is aircraft mass, Q is dynamic pressure, and S is wing area. The 
thrust is determined from tabular data that relate actual thrust to the particular 
engine parameter recorded. A second set of expressions for the force coefficients 
is obtained by specifying the corresponding aerodynamic models of the form 


C L = C L (a> M) + £*i (7.4) 

i 


C D = C D (a ’ M) + ^ d i (7.5) 


5b 


(7.6) 


c c = c c (s, M) + £c. 

where M is the Mach number. Terms forming the sums in equations (7.4) through 
(7.6) represent the contributions of angular rates, flaps, spoilers, control sur- 
faces, landing gear, and ground effects. It should be emphasized that, in general, 
the coefficient models represent all that is known about the aerodynamic properties 
of the aircraft from theoretical predictions, wind-tunnel experiments, and flight 
testing. 

To estimate an angle-of-attack time history, the lift-coefficient expressions 
of equations (7.1) and (7.4) are equated, giving a nonlinear algebraic equation to 
be solved for angle of attack at each time point. An iterative procedure like the 
Newton-Raphson method works very well for this problem. A similar technique is used 
with side-force coefficient expressions for estimating the sideslip angle. The 
coefficient method yields good, "wide-band" estimates of both angles of attack and 
sideslip (ref. 53). In the analysis of the Flight 191 data described later in the 
chapter, however, measured flow angles (left and right alpha vanes) were used to 
derive angle of attack, whereas the coefficient method was used to estimate the 
angle of sideslip (beta-vane measurements are not included with DFDR records). 


Wind Estimation 

As illustrated in chapter 4, the SMACK procedure can be used with data from 
several sources (e.g., flight recorder and ATC radar) to determine the wind pattern 
along the flightpath of an aircraft. This technique has been useful in the analysis 
of recent airliner encounters with severe turbulence. To solve the aircraft flight- 
path wind problem discussed in the next section, one procedure would use SMACK to 
determine the integrator initial conditions, accelerometer biases, and forcing- 
function time histories t d x , dy, d^), (<J>, 9, >l>), and ( g x , g„, g h ) that provide the 
"best fits" to the measurement records (x, y, h), (<t>, 0, 4>), (V, a, 6), and (a x , ay, 
a ) (i.e., those that minimize eq. (2.3)). The wind estimates (W X y, V wc j) along 

the flightpath are obtained as part of the SMACK solution. 

A second procedure would fit only the inertial data using SMACK and then calcu- 
late the wind components separately from 

w = x - V cose cosil) 
x w w 

w y = y - V cose w siml> w 
w h = h - V sin9 w 

where the wind-axis Euler angles (9 W , l w ) are given by 

6 = sin"^(cosct cose sin0 - C cose) (7.10) 

w 


(7.7) 

(7.8) 

(7.9) 
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(7.11) 


* 


w 


i|» + tan 


-1 


sing cos4) - sina cose sin4>)/D] 


where 

C = sina cose cos4> + sine sin<(> ; D = cosa cose cos0 + C sine 

Because the air-data variables (V, a, 6) derived from the Flight 191 DFDR were 
relatively "smooth," the second procedure was used to obtain the wind estimates 
presented in the next section. The coding list for this application is shown in 
figure 7.3. 


The Flight 191 Accident 


The Delta Airlines Flight 191 windshear accident occurred during an attempted 
landing at the Dallas/Ft. Worth Airport. A 14-min portion of the approach path of 
the L— 1011 aircraft as measured by ATC radar is shown in figure 7.4(a). The figure 
shows the aircraft approaching from the northeast and turning south onto the final 
glidepath to the runway. The point of initial contact with the ground (shown as a 
solid circle in fig. 7.4) is at a location 361 ft to the east and 6343 ft to the 
north of the runway origin. Figure 7.4(b) shows for the same time interval, two 
curves of pressure altitude that represent the "best fu" of the DFDR record to the 
transponded record. The time shift required to synchronize the two data sources was 
one second. 
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Figure 7.3.- Coding i . st. for the SMACK state-estimation procedure. 
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Figure 7.4.- Flight 191 approach, (a) ATC radar groundtrack, (b) transponded and 
DFDR altitude records compared for time synchronization. 

Following data synchronization, and calculation of true airspeed, angle of 
attack (from the vanes), and sideslip angle (from the performance equations), the 
state-estimation program SMACK was applied to determine the winds. The kinematic 
equations were integrated over a five-minute period that starts before the turn onto 
final approach and ends with the initial ground contact. The fits to the position 
data for this period are shown in figures 7.5(a) and (b). In these figures the 


58 




small circles represent the measured values and the dashed lines represent the time 
histones generated by SMACK. When the fact that the least count of the ATC 
tracking data is on the order of 1/8 nautical mile is taken into consideration, 
there is good agreement between the estimated path and the radar groundtrack in 
figure 7.5(a). The inertial altitude estimate is compared with the DFDR barometric 
altitude record m figure 7.5(b). During most of this five-minute interval, there 
is good agreement between the estimated inertial altitude and the measured 
barometric altitude. However, during the final portion (in the downburst) there is 
some discrepancy, apparently due to local pressure variations caused by the 
atmospheric disturbance. 


Because the aircraft was in the downburst for less than a minute before its 
initial contact with the ground, the rest of the analysis described in this section 
will cover only the final 60 seconds of flight. Figure 7.6 shows time histories of 
the three body-axis accelerations, while figure 7.7 shows time histories of the 
three body-axis Euler angles. Mote that the plots of figures 7.6 and 7.7 include 
the SMACK-derived "best fits" to the DFDR data records. Figure 7.8 shows time 
istories of the aerodynamic variables (true airspeed, angles of attack and side- 
s ip). The angle of attack was computed from the average value of right and left 
vanes after correction for upwash and pitch rate. Since vane-rate limiting (at 
about 19°/sec) occurred during the last 20 sec, the rapid excursions in angle of 
attack shown in figure 7.8(b) are probably attenuated. As mentioned earlier the 
angle of sideslip was computed from the measured side force using predicted aerody- 
namics and including terms for rudder deflection and yaw rate. 

Figure 7.9(a) presents a time history of the aircraft heading angle shown with 
the groundtrack angle. The value observed for the groundtrack angle at the final 
ime is 17 0 from true north. This estimate of groundtrack angle at the final time 
is m agreement with the orientation of the landing gear marks found in the field 
where the aircraft first contacted the ground. Figure 7.9(b) presents a time his- 
tory of the true airspeed together with the estimated groundspeed. The groundspeed 
is seen to be increasing beyond 210 knots at the point of initial contact. Fig- 

U k 6 I'L 3 !! w S ^ during the final few seconds there appears to be a tailwind of 
about 60 ft/sec (35 knots). 

.. ,^ he general Pattern of the winds can be deduced from figure 7.10, which shows 
the three components of the wind vector. The horizontal components are shown in 
figure 7.10(a); the vertical component is shown in figure 7.10(b). Because the 
vertical wind estimate depends on the angle of attack, the vane rate-limiting men- 
loned earlier will also attenuate the vertical wind excursions. The results shown 
in figure 7.10 indicate that the aircraft encountered a strong downflow for a time 
period of 20 sec followed by a rapid change in vertical wind direction, followed by 
further changes about 5 sec apart. During the period of major downflow, the air- 
craft experienced vertical winds on the order of -10 to -40 ft/sec. When the air- 
" a f , t , entere J, the downflow, the headwind increased from about 20 ft/sec to more than 
50 ft/sec. Then, during a period of 26 sec, there was a change to a tailwind of 
more than 50 ft/sec. ° 
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Figure 7.5.- Position solutions (5 min), (a) Groundtrack, (b) pressure altitude. 

Figure 7.11 shows winds along the flightpath from different perspectives that 
clearly indicate the pattern of winds in the downburst. The diagram in fig- 
ure 7.11(a) shows the flightpath viewed from above with the wind arrows computed 
from the horizontal components w x and w . These results show the changes in the 
magnitude and direction of the horizontal wind as the aircraft proceeds through the 
downburst. As shown by the rotation of the horizontal wind vector, the source of 
the downflow appears to be located west of the flightpath. The diagram in 
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Figure 7.9.- Velocity estimates (last 60 sec), (a) True heading and groundtrack 

angles, (b) true airspeed and groundspeed. 

figure 7.11(b) shows the flightpath viewed from the west with the wind arrows 
computed from the w x and w h components. Following the downflow portion, the 
outflow near the ground is evident along with changes in the vertical wind. The 
before and after the downflow indicate the presence of vortex rings 
(ref. 54). According to the vortex-ring model, when a ring impacts the surface, its 
circulation is spun up, providing a mechanism for the changes in vertical wind that 
are observed near the ground. In particular, the rapid changes just after the 
downflow are typical of a series of strong vortices. 

The analysis of this accident, which is one of the first involving an aircraft 
with an onboard digital recording system, provides a detailed look at the pattern of 
low-level windshear in a downburst environment. The winds derived from the Delta 
191 records provide important new information to augment ongoing experiments and 
theoretical research on the downburst phenomenon. Furthermore, the estimated wind 
time-histories provide a new set of data that represents a severe downburst for use 
in simulation and pilot training. 
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Figure 7.10.- Flightpath wind-component estimates (last 60 sec), (a) Southerly and 

westerly, (b) vertical. 
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Figure 7.11- Flightpath wind vectors (last 44 sec), (a) Seen from above 

(b) seen from the west. 



8. CONCLUDING REMARKS 


In this report the evolution of state-estimation technology for use in the 
analysis of aircraft flight data has been traced. The underlying mathematics have 
been reviewed, and a general-purpose aircraft state-estimation program called SMACK 
as been described. Three recent applications involving the estimation of winds 
uler angles, and air variables have been discussed, and several examples based on 
the applications have been presented. The examples demonstrate that a general- 
purpose state estimation program can be used to solve the applications discussed, 
n one additional example of an application not previously reported, it was shown 
that inertial position measurements and onboard strap-down measurements can be 

combined using a state-estimation procedure to provide estimates of Euler angles and 
air-data variables. 


The coding procedure .or solving flight-data consistency problems with SMACK 
has been introduced, and two "real" flight-data applications have been discussed 

G JT’ f u fll S ht - test methodology for identification of an aerodynamic model 
that includes the use of state estimation was described. In the second, the appli- 
cation of state estimation in the analysis of a windshear accident was presented. 

is hoped that the text portion of this report has helped to make the flight-data 
analyst aware of the potential advantages of using state estimation in solving a 
variety of problems, and that it and the appendices that follow will serve ade- 
quately as a User's Manual for the SMACK computer program. 
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APPENDIX A 


PROGRAM DESCRIPTION 


This appendix presents an overview of SMACK which is intended to aid the user 

m R ™T^ and , ing Pr ° gram Strucfcure and flow - T he program is written entirely in 
FORTRAN 77 and consists of approximately 8000 lines of code (5600 statements) The 
main program and two of the most important subroutines will be described in some 
detail here. However, every routine used by SMACK will at the least be mentioned 
and its place m the calling hierarchy will be described. For ease of reference ’ 
the calling hierarchy is shown in figure A1. Further information on program imple- 
mentation will be found in appendix B; output listings for test problems are 

mrA 10 appandlces C and D ’ and instructions for preparing the user subroutine 
DATA are presented m appendix E. 


The Main Program (SMACK) 

A block diagram of the main program is shown in figure A2. It calls several 
subroutines that initialize arrays, read and analyze the problem coding list 

rr rd ’ Pr ° Vide a startin g trajectory, solve the problem, and ’print 
and plot the solution. The subroutines called directly by SMACK are described as 


INIT initializes all arrays in the COMMON blocks. Maximum dimensions and 
values for constants are set in a BLOCK DATA routine. 

READ accepts and analyzes the problem coding list, sets solution, estima- 
tion, and measurement flags, and determines entries in the Jacobian 
matrices. The "bookkeeping" chores performed here allow efficient manip- 
ulation of sparse matrices later in the program. 

ARIN displays the contents of state-variable, forcing-function, and 
output-variable integer arrays, including the Jacobian matrix struc- 
ure. The information may be helpful for understanding how SMACK works. 

M0DL creates an aircraft flight trajectory consisting of a rising 
coordinated 180° turn in wind, adds noise (if specified) to the measure- 
ment records, and stores these records, along with the "true" time his- 
tories of all output variables for later comparison with estimates. The 
trajectory is used for program testing. 

DATA accesses an external flight-data record for analysis by SMACK. The 
preparation of this subroutine, which must be supplied by the user is 
discussed in appendix E. ’ 
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COMMENTS 


CALL IIM IT 
CALL READ 



INITIALIZE ARRAYS, READ AND ANALYZE 
CODING LIST, SET UP PROBLEM 

OUTPUT SUMMARY OF CODING LIST 


OUTPUT CONTENTS OF INTEGER ARRAYS 


OBTAIN DATA RECORD FROM DATA 
(USER-SUPPLIED), OR MODL (A/C 
SIMULATION) 

OUTPUT SUMMARY OF SIMULATED 
MANEUVER 

OBTAIN INITIAL CONDITIONS, FORCING- 
FUNCTION TIME HISTORIES NEEDED TO 
GENERATE THE STARTING TRAJECTORY 

OUTPUT SUMMARY OF STARTING 
SOLUTION 


OUTPUT PLOTS OF STARTING SOLUTION 


OUTPUT CONTENTS OF TIME ARRAYS AT 
t = 0, t = DT 


Figure A2.- Flow diagram of SMACK (main program). 
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OBTAIN AN IMPROVED STARTING 
SOLUTION, EVALUATING PARTIAL 
DERIVATIVES FROM MEASUREMENTS 



OUTPUT SUMMARY OF STARTING 
SOLUTION 

OBTAIN NORMAL SMACK SOLUTION 
OUTPUT SUMMARY OF FINAL SOLUTION 


PERFORM ACCIDENT-ANALYSIS CHORES 

OUTPUT SUMMARY OF ACCIDENT 
ANALYSIS 


OUTPUT PLOTS OF FINAL SOLUTION 


RETURN TO DATA (FOR USER), OR TO 
MODL FOR SOLUTION TABULATION 

OUTPUT SIMULATION RESULTS 


Figure A2.- Concluded. 
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STRT derives the initial conditions and forcing-function time histories 
necessary to create the starting nominal trajectory required by the SMACK 
algorithm. It also determines the diagonal elements of the weighting 
matrices used in the performance measure. 

VRWC calculates and displays mean and rms values of measurement-residual 
and forcing-function time histories after a solution has been obtained. 

PRNT displays values of initial conditions, bias errors and scale factors 
after a solution has been obtained. This routine is also called by MODL. 

ARNM displays the contents of state-variable, forcing-function, and 
output-variable real arrays, as well as the Jacobian matrices at che 
first time point calculated. This information may help in problem 
diagnosis. 

SECOND (a Cray utility) determines the CPU time for the computationally- 
intensive parts of the program. It is included with the VAX versions to 
work with the LIB$INIT TIMER VAX utility. 

MINC performs the minimization of the cost function (performance mea- 
sure), utilizing the backward-filter, forward-smoother algorithm dis- 
cussed in chapter 2, and outlined in table 2.1. 

SOLU completes accident analyses by computing Euler angles from position 
data (when appropriate) and prints a summary of results, which includes 
groundspeed, groundtrack, airspeed, flightpath angle, angle of attack 
magnetic heading, lift and drag forces. 

PLOT generates printer or DISSPLA plots of output-variable time 
histories . 


he several places at which output is produced during problem solution are 
indicated in the block diagram of figure A2. Notice that results will be printed 
for the starting solution, for an improved starting solution (if selected), and for 
the final solution. Plots can be made of either the starting solution or the final 
S r i 1Cm ‘ SOlution °P tion ^ags ISTRT , ISOLU, IDATA, IPLOT, and IDBUG are spec- 

lf ; e the probiem codln g list covered in chapter 5. The first three of these are 
set by the solution description statement (ENG or MKS), with 


ISTRT = J, ISOLU = K, IDATA = L, 

while flags IPLOT and IDBUG are set by the END statement, with 
IPLOT = I, IDBUG = J. 
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The Starting Routine ( STRT ) 


The SMACK program requires a starting trajectory "close enough" to the optimum 
oath to provide reasonable assurance that the solution algorithm will converge The 
initial conditions and forcing functions necessary to create a starting trajectory 
are derived by subroutine STRT. A block diagram for carrying out these important 
calculations is shown in figure A3. The procedure begins with the filtering of each 
measurement record, which provides estimates of residual variances for use in t e 
measurement-error weighting matrix R in equation (2.3). Algebraic methods are 
then used to determine state-var iable estimates from the available measurements, 
from which are obtained the required initial conditions and forcing functions an 
corresponding variances for elements of the weighting matrix Q in equa- 
tion (2.3). Subroutines called directly by STRT are the following. 


ATRK eliminates all 360° jumps in heading records, thus creating 
continuous-angle time histories. This routine is also called by TRAJ and 

ACSIM. 


FILT provides digital filtering of a record with estimates of the first 
and second time derivatives. The filter frequency characteristic is that 
of a fourth-order, zero phase-shift, low-pass filter (see appendix ). 

AVAR computes the sample mean and variance for a time history. This 
routine is also called by MODL and VRNC. 


COOR computes position coordinates from range, bearing, and elevation 
data and the location of the tracking site. 

CINS computes vehicle inertial velocity components from groundspeed and 
groundtrack data* 

FOIL computes roll and pitch angles from airspeed, heading, altitude, 
winds, and aircraft performance data (ref. 47). 

RADR computes roll, pitch, and yaw angles from radar position winds, and 
aircraft performance data (ref. 47). This routine is also called by 
ACSIM and SOLU. 


ANGL computes roll, pitch, and yaw angle time derivatives from rate-gyro 
data and integrates them to obtain the Euler-angle time histories. 

WIND computes wind components in the Earth frame from wind magnitude and 
direction data. 


AERO computes components of vehicle velocity in the Earth frame, with 
respect to the air mass, from air and Euler-angle data. 


78 



ENTER 


COMMENTS 



THE FILTERING ROUTINE 
FITS THE RECORD AND 
ESTIMATES THE FIRST 
TWO DERIVATIVES (SEE 
APPENDIX E) 


COORDINATES GIVEN BY 
x = x f + R cos B cos E 
y = Y r + R sin B cos E 
h = h r + R cos B sin E 


EULER-ANGLE ESTIMATES 
OBTAINED FROM AIRSPEED 
HEADING, ALTITUDE, OR, 
FROM POSITION DATA [47] 


ANGLES DETERMINED BY 
INTEGRATING ANGULAR 
RATES FROM 

4/ = (qsin <j> + rcos </>)/ cos 6 
6 = qcos0 - rsin <p 
<p = p + ip sin 8 


Figure A3.- Flow diagram for subroutine STRT. 
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WINDS 



WIND COMPONENTS ARE 

«x ■ - w xy cos W hd 
w y = - w xy sin W hd 
w h ~~ ^wd 


VELOCITY COMPONENTS ARE 
u = V cos a cos j3 
v = V sin 0 
w = V sin a cos|3 


K1 


u 

Va 

= L t 

V 

Ua. 


w 


* = *a + w x 

y = y a + w v 
h = h a + w h 


VELOCITIES DETERMINED 
BY INTEGRATING INERTIAL 
ACCELERATIONS FROM 


x: 
i 




o 

1 

Y 

= L* 

a y 


0 

h 




CQ 


Figure A3- - Concluded. 
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VELS computes components of acceleration in the Earth frame from angle and 
accelerometer data and integrates them to obtain inertial-velocity time 
histories. 


The Minimizing Routine (MINC) 

The subroutine that directs the procedure of minimizing the performance measure 
is MINC. The function of this subroutine is to carry out the steps of the SMACK 
algorithm outlined in table 2.1. Its block diagram is shown in figure A4, where it 
can be seen that when the performance measure (cost) is reduced (or first computed), 
changes in initial conditions and forcing— function time histories are calculated, 
and the iteration counter is tested. If the iteration count is less than NIT, an 
update occurs. If the cost is not reduced by this update, the change is halved and 
the update is recalculated. If this process succeeds in reducing cost, another 
iteration is performed; if not, it terminates with a message "NO IMPROVEMENT ON THIS 
ITERATION". The subroutines called by MINC are described as follows: 

TRAJ solves the differential equations of the state model and evaluates 
the measurement model (see fig. 3.1), and the performance measure. Note 
that this is the first step of the SMACK algorithm outlined in table 2.1. 

SENS solves the backward-information filter, calculates the parameter 
changes, and solves the forward smoother to determine the forcing- 
function changes, which are the steps of table 2.1(b). When the solution 
is complete (IFIN=1), it calculates estimates of parameter standard 
deviations from the diagonal elements of the information-matrix inverse. 

SETP performs the parameter and forcing-function update (the last step 
outlined in table 2.1). 


Other Subroutines 

For the rest of the discussion concerning program structure, it will be con- 
venient to refer to the subroutine calling hierarchy shown in figure A1 . The sub- 
routines yet to be mentioned are first called by READ, MODL, ARNM, SOLU, and PLOT 
and will be considered in that order. They are as follows: 

JACSX sets up the three Jacobian (partial derivative) arrays FX, FW, and 
HX according to the problem to be solved. 

ACSIM computes the required initial conditions and forcing-function time 
histories for an aircraft maneuver consisting of a rising, coordinated 
180° turn in a wind environment. 
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EVALF evaluates the state-variable time derivatives (see the the state 
model of eq. (2.1) and fig. 3.1). This routine is called by MODL, ARNM 

anH TRfl I J ’ * 


EVALH evaluates the output variables (see the measurement model of 
eq. (2.2) and fig. 3.1). This routine is called by MODL, ARNM, and TRAJ 

NOYZ creates a sequence of Gaussian random numbers with user-specified 
variance. The sequence is used in MODL to contaminate simulated measure 


GGNQF (An IMSL routine) chooses a number at random from a "Gaussian 
hat" — a pseudo-random normal (0, 1) deviate. 

TRANS transforms user-specified measurements for use as forcing func- 
tions, under an option covered in chapter 5 that is useful in some data- 
consistency experiments: rate-gyro measurements are transformed to 

Euler-angle rates (eq. (3.12)); accelerometer measurements are trans- 
formed to accelerations in the Earth frame (eq. (3.11)); wind measure- 
ments are transformed to components in the Earth frame. This routine is 
called by ARNM and TRAJ. 


EVLFX evaluates active elements of the Jacobian array FX (partial- 
is™^ matriX defined in eq - (2 - 6)) - This r ° utine is called by ARNM 


EVLFW evaluates active elements of the Jacobian array FW (partial- 

ann^ ” atrlX defl " ed ^ 6q ' <2 ' 6)K ™ S r ° Utine 13 Called b * 


EVLHX evaluates active elements of the Jacobian array HX (partial- 
is™^ 6 matriX defin6d in ec ^ (2 - 6 >-> This routine is called by ARNM 


LEQT2P (an IMSL routine) computes the solution of the linear equations 
imp led by equation (2.13) to obtain a parameter change DELP. The infor- 
mation matrix EM is stored in symmetric-storage mode. This routine is 
called each iteration by SENS, when IFIN=0 


NV3P (an IMSL routine) computes the inverse of the information matrix 
to obtain estimates of parameter standard deviations. This routine is 
called by SENS when a solution is complete (IFIN=1). 


CONV tr^sforms accelerations from body to wind axes to make calculations 
lift and drag forces; computes groundspeed, groundtrack, flightpath 
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angle, indicated airspeed, and magnetic heading (all variables of inter- 
est for accident analysis). 

PGPLO produces a printer plot with the time axis running the length of a 
page. A total of 300 points per variable may be displayed. 

USPLO (an IMSL routine) produces an x-y printer plot, one plot per 
page. 

XYPLO produces x-y plot files with DISSPLA for output by an external 
device (installation-dependent). 
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APPENDIX B 


PROGRAM IMPLEMENTATION 


This appendix considers some aspects important to program implementation, such 
as the computer systems utilized, the structure of COMMON and the way in which 
calculations involving the Jacobian matrices are performed. Should the user desire 
to make modifications to the program, an example showing how to add a new output 
variable is included. Finally, some options are suggested for providing the large 
temporary storage required in subroutine SENS, since such storage may be 
machine-dependent. 


Computer Systems 

SMACK has been implemented at Ames Research Center on the Cray X-MP, Y-MP, and 
VAX 11/785 and 8650 computers. Running times on the VAX 8650 (double-precision) 
version have been observed to be about an order of magnitude greater than the Cray 
for the same problems. The user should note that the Cray word is 64 bits (8 
bytes), while the single-precision VAX word is 32 bits (4 bytes). Differences 
between the Cray and VAX versions are few, and are indicated in the code with the 
comments "CRAY SPECIFIC" or "VAX SPECIFIC" and "VAX DP". The double-precision VAX 
version was created by adding to the beginning of each routine the statement 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 'VAX DP 
and accessing the double-precision IMSL library. 

The program has been designed to be portable. For example, all Hollerith 
variables are stored four characters per word, a feature compatible with most 
machines. The SMACK program layout is shown in figure B1, and consists of the 
COMMON groups, the main program, and the subroutines described in appendix A Those 
routines that make use of the IMSL or DISSPLA libraries are indicated. The terms 
*C0MDECK and "DECK are directives to the UPDATE utility, the Cray program "librar- 
ian". The directives are considered as comments in the VAX versions of SMACK. It 
should be noted that the VAX double-precision version has been submitted to 
COSMIC. Computer Software Management and Information Center, Suite 112, Barrow 
Hall, The University of Georgia, Athens, GA 30601 (404-542-3265). 


COMMON Structure 

The COMMON area used by SMACK is divided into three groups, C0MTME , COMVAR and 
COMNME . The first group, C0MTME , shown in table B1, includes all the time-history 
arrays for the program, and requires 1,350,000 words of memory to store 6000-point 
records. In block /STM/, array W will contain forcing-function time histories, from 
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FILE STRUCTURE 


COMMON DECKS: *COMDECK COMTME 

★COMDECK COMVAR 
*COMDECK COMNME 


MAIN PROGRAM: 


*DECK SMACK: 

SMACK 


SUBROUTINES: 



*DECK SETUP: 

INIT r 

READ, JACSX 

*DECK START: 

STRT , 
COOR, 

AERO, ANGL, ATRK, CINS, 
FOIL, RADR, VELS, WIND 

★DECK OPTIM: 

MINC , 

SENS, SETP , TRAJ 

*DECK EVLMD : 

EVALF , 

EVALH, TRANS 

*DECK EVLJC: 

EVLFW, 

EVLFX, EVLHX 

★DECK ACMOD : 

MODL, 

ACSIM, NOYZ 

★DECK DEBUG: 

ARIN, 

ARNM 

★DECK MISCL : 

AVAR, 

CONV, FILT, SECOND (1) 

★DECK PRNTS : 

PRNT, 

SOLU, VRNC 

★DECK PLOTS: 

PLOT , 

PGPLO, XYPLO (2) 

EXTERNALS : 



IMSL LIB (3) : 

GGNQF 

, LEQT2P, LINV3P, USPLO 

USER CREATED: 

DATA 



PROGRAM NOTES 

(1) SR SECOND A CRAY UTILITY (SUPPLIED FOR VAX VERSION 
TO BE USED WITH VMS LIB$ INIT_TIMER UTILITY) . 

(2) SR XYPLO UTILIZES THE DISSPLA LIBRARY FROM ISSCO, 
10505 SORRENTO VALLEY ROAD, SAN DIEGO, CA 92121. 

(3) IMSL INC., 7500 BELLAIRE BLVD . , HOUSTON, TX 77036. 


Figure B1.- SMACK file structure and program notes. 
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TABLE B 1 . - *COMDECK COMTME 


COMMON/STM/X(21 ,6000) ,W(9,6000) ,DMW(54000) 
DIMENSION DW(9,6000) 

EQUIVALENCE (DW,DMW) 

DIMENSION XA(6000,3) , YA(6000,3) 

EQUIVALENCE (XA,DMW( 1 ) ) , (YA,DMW( 18001)) 
COMMON/MSM/Z(30,6000) ,V(30,6000) ,KZN(30,6000) , 

* H(30, 6000) ,D(30, 6000) , S(30,6000) 
DIMENSION Y( 30, 6000) 

EQUIVALENCE (Y,V) 

COMMON/ AUX/XN( 6000) , YN( 6000) , ZN( 6000) , FN(6000) , 

* KFN(6000),NFN(6000) 


the list of table 5.2, while X will contain state-variable time histories from the 
list of table 5.3. The array DW contains the forcing-function updates, while XA and 
YA are used for plotting. In block /MSM/ , arrays Z and KZN are used for storage of 
measurement time histories and H for output-variable time histories, in the order 
specified by table 5.1 (see also appendix C). The array V contains residuals for 
performance measure evaluation, and arrays D and S are used in STRT to store 
derivative estimates, and again later to store other time histories needed in the 
SMACK solution. The use of array Y is discussed in appendix E. Included in block 
/AUX/ are six 6000-word scratch arrays. 

The second group, C0MVAR, is shown in table B2. It defines the rest of COMMON 
memory, and requires about 5,760 words. The equivalencing for variable array ED 
indicated by table B2 is illustrated in figure B2. The array ED has component 
arrays WD, XD, and HD which are used to hold single time-point values of forcing 
function W, state X, and output H arrays, respectively. The overlap among the com- 
ponent arrays occurs because six variables (PH2, TH2, PS2) and (X2, Y2, H2) can be 
either forcing functions or state variables, while another six variables (PHI, THT, 
PSI) and (X, Y, H) are both state variables and output variables. A similar equiva- 
lencing exists for the name array NME and its components (NMW, NMX, NMH), the units 
array UNE and (UNW, UNX, UNH), the scaling array SCE and (SCW, SCX, SCH) and the 
index array KE and (KW, KX, KH). 

The equivalencing scheme employed for the parameter vector XP in table B2 is 
also shown in figure B2. Array XP has component arrays X0, WB, VB, and SF which 
represent the initial conditions (table 5.3), forcing-function means (table 5.2) 
bias errors, and scale factors (table 5.1), respectively. Although the length of 
the parameter vector XP is 96, the maximum number of parameters that can be iden- 
tified is 42. The parameter addresses are stored in array IXP. Notice that there 
is no overlap among the parameter components. A similar equivalencing exists for 
the "statistics" vector SXP with components (SX0, SWB, SVB, SSF), and the parameter 
index vector KXP and its components (KX0, KWB, KVB, KSF). 

The third COMMON group C0MNME is shown in table B3. This group consists 
entirely of EQUIVALENCE statements which make possible the use of variable names in 
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TABLE B2 . - *COMDECK COMVAR 


COMMON/MOD/ ALFO , CLAO , WGLD , V ARM 
COMMON/ DEF/NTITL( 12) ,SKL(8) ,UNM( 12) , NIT 
COMMON/TIM/DT , NPTS , NINT , NSMP , NSKP , MPTS , FC 1 
COMMON/SKL/CMK , CMF , CRD , CNL , CMM , CMG 
COMMON /CON / RG AS , TSSL , RHOO , P I , C ( 30 ) 

COMMON/DEC/ IH1 (6) , IH2(6) , IS0(6) , IS1 (6) , IS2(6) ,110(6) ,111(6) 
COMMON/FFN/IW( 15) ,NW, IWU( 15) ,NWU,NUT,W0(9) ,DW0(9) ,NWT, 

* KIN(15),HN( 15) ,NIN,WU( 15) ,WP( 15) 

C0MM0N/STV/IX(21 ) ,NX,NXT,KSX(21 ) , ISX(21 ) ,NSX, 

* KF(21),IF(21),NF,FD(21), 

* KZX(21 ),IZX(21),NZX,IXZ(21) ,NXZ 
COMMON /OUT / KZ ( 30 ) , IZ( 30) ,NZ , IH( 30) ,NH,NHT,NST, 

* KV(30),IV(30),NV,VD(30),DD(30),SD(30), 

* IZH(30) ,NZH, I HZ (30) ,NHZ, IZS( 15) ,NZS, 

* HO(30),DO(30),SO(30),PLS(30,5),FC(30),NZPT(30) 
COMMON/ JAC/LFX( 72),MFX( 72),KFX( 72),FX( 72) ,NFX,NFXT, 

* LFW( 15) ,MFW( 15) ,KFW( 15),FW( 15) ,NFW,NFWT , 

* LHX( 196) , MHX( 196) ,KHX( 196),HX(196),NHX,NHXT, 

* LHTH( 196) ,MHTH( 1089) ,KPHX( 1089) , 

* MPFW( 42 , 15 ) MPFX(42,72) 

COMMON /VNC/WS( 15) ,WMN( 15) ,WSD( 15) ,QD( 15) ,QI( 15) , 

* VS(30),VMN(30),VSD(30),RD(30),RI(30) 
COMMON/PAR/KXP(96),IXP(42),NXP,NXPT,XP(96), 

* SXP(96) ,XPB(96) ,POI(96) ,DELP(42) ,NPAR, 

* I WB ( 1 5 ) , NWB , NWBT , I VB ( 30 ) , NVB , NVBT , ISF ( 30 ) , NSF , NSFT 
COMMON/ FLG/KANG , KXYH , KR AD , KWND , K AI R , KSPF , KPQR , KMOM , KPOS , KVEL , 

* IANG, IXYH, I RAD, IWND, I AIR, ISPF, IPQR, IMOM , IPOS, IVEL, 

* LWND, LAIR, LSPF,LPQR,LMOM,IAC, IPS, IVN, IRE, IWD, 

* I SOLU , I DAT A , I DBUG , ISTRT , I PLOT , INOYZ 
COMMON/MST/KE ( 54 ) , I E ( 54 ) , NET , NME ( 54 ) , ED( 54 ) , 

* SCE(54),ISCE(54),UNE(54),IUNE(54) 

DIMENSION WD( 15), KW (15) ,SCW( 15) ,UNW( 15) ,NMW( 15) , 

* WB( 15) ,KWB( 15) ,SWB( 15) , 

* XD( 2 1 ) , KX( 2 1 ) ,SCX(21 ) ,UNX(21 ) ,NMX(21 ) , 

* X0(21 ) ,KX0(21 ) ,SX0(21 ) , 

* HD(30) ,KH(30) ,SCH(30) ,UNH(30) ,NMH(30) , 

* VB(30) ,KVB(30) , SVB ( 30 ) , SF ( 30 ) , KSF ( 30 ) , SSF ( 30 ) 
EQUIVALENCE ( WD, ED(1)),( XD, ED(10)),( HD, ED(25)), 

( KW, KE ( 1 ) ) , ( KX, KE( 10) ) , ( KH, KE(25)), 
(SCW,SCE( 1 ) ) , (SCX,SCE( 10) ) , (SCH,SCE(25) ) , 
(UNW,UNE( 1 ) ) , (UNX,UNE( 10) ) , (UNH,UNE(25 ) ) , 
(NMW,NME( 1 ) ) , (NMX ,NME( 10) ) , (NMH,NME(25) ) 
EQUIVALENCE ( XO, XP(1)),( WB, XP(22)),( VB, XP(37)), 

* ( SF, XP(67 ) ) , 

* (KXO , KXP( 1 ) ) , ( KWB ,KXP(22 ) ) , (KVB,KXP(37) ) , 

* ( KSF , KXP( 67 ) ) , , „ 

* ( SXO , SXP ( 1 ) ) , ( SWB , SXP ( 22 ) ) , ( SVB , SXP ( 37 ) ) , 

* (SSF,SXP(67 ) ) 
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TABLE B3.- *COMDECK COMNME 


EQUIVALENCE (DL ,ED( 1)) 
(DX , ED( 4)) 
(GX ,ED( 7)) 

( PH2 , ED( 10) ) 
(X2 ,ED( 13) ) 
(WX , ED( 1 6 ) ) 
(PHI ,ED( 19) ) 
(XI ,ED(22) ) 

( PHI , ED( 25 ) ) 
(X ,ED(28) ) 
(RNG,ED(3D) 
(WXY,ED(34)) 
(VT ,ED(37) ) 
(AX ,ED(40) ) 
(P ,ED(43) ) 
(AL ,ED(46) ) 
(RNA,ED(49) ) 
(HDG,ED(52) ) 
EQUIVALENCE (PH2D,FD( 1) 
( X2DT , FD( 4) 
(WXDT,FD( 7) 
( PHI D , FD( 10) 
(X1DT,FD( 13) 
(PHID,FD( 16) 
(XDOT,FD( 19) 
EQUIVALENCE (UA ,SD( 1) 
( AXB , SD( 4) 
( UAPS , SD( 7) 
(UAVN,SD( 10) 
EQUIVALENCE ( X 1 A , DD( 10) ) 
(VT1 ,DD( 13) ) 
EQUIVALENCE (XAC,C( 1)), 
(XPS,C( 4)), 
( XVN , C( 7)), 
(RX ,C(10)), 
> (RXA,C( 13) ) , 


, (DM ,ED( 2 ) ) , ( DN ,ED( 3)), 
,(DY ,ED( 5 ) ) , ( DH ,ED( 6)), 

, (GY , ED( 8 ) ) , ( GH ,ED( 9)), 

, (TH2 ,ED( 1 1 ) ) , ( PS2 , ED( 12 ) ) , 

, ( Y2 ,ED(14)),(H2 ,ED(15)), 
,(WY ,ED( 17) ) i (WH ,ED( 18) ) , 

, (TH1 ,ED(20) ) , ( PS1 , ED(21 ) ) , 

, ( Y1 ,ED(23)),(H1 ,ED(24) ) , 

, (THT ,ED( 26 ) ) , ( PSI , ED(27 ) ) , 
,(Y ,ED(29)),(H , ED(30) ) , 

,(BRG,ED(32)) t (ELV,ED(33)), 
,(WHD,ED(35)),(VWD,ED(36)), 
,(AV ,ED(38)),(BV , ED(39) ) , 

, ( AY ,ED(4l ) ) , ( AZ ,ED(42) ) , 

, (Q ,ED(44)),(R ,ED(45) ) , 

, ( AM ,ED(47) ) , (AN ,ED(48)) f 
, (BRA, ED (50) ) , (ELA,ED(51 ) ) , 

, (VGR, ED(53) ) , (TRK,ED(54 ) ) 

) , (TH2D , FD( 2) ) , (PS2D,FD( 3)) 
),(Y2DT,FD( 5) ) , (H2DT,FD( 6)) 
) , (WYDT,FD( 8) ) , ( WHDT,FD( 9)) 
),(TH1D,FD(11)),(PS1D,FD(12)) 
),(Y1DT,FD(14)),(H1DT,FD(15)) 
),(THTD,FD( 17)),(PSID,FD(18)) 
) ,(YD0T,FD(20)),(HD0T,FD(21)) 
),(V A ,SD( 2) ) , (WA ,SD( 3)) 
) , ( AYB ,SD( 5) ) , ( AZB ,SD( 6)) 
) , (VAPS,SD( 8) ) , (WAPS,SD( 9)) 
) , ( VAVN,SD( 11)), (WAVN,SD( 12) ) 
,(Y1A,DD(11)),(H1A,DD(12)), 

, (AVI ,DD( 14) ) , (BV1 ,DD( 15) ) 

( YAC , C( 2) ) , (ZAC,C( 3)), 
(YPS,C( 5)),(ZPS,C( 6)), 
(YVN,C( 8) ) , (ZVN,C( 9)), 

(RY ,C( 1 1 ) ) , ( RH ,C( 12) ) , 
(RYA,C( 14) ) , (RHA,C( 15) ) 


subroutines that carry out specific model-related calculations (such as in EVALF and 
EVALH). Note that the array SD holds single time-point values of time histories 
stored in the S array. An alphabetical listing of all SMACK routines summarizing 
their COMMON requirements is given in table B4. Notice that each routine that 
requires COMTME or COMNME also requires COMVAR, but that no routine uses both COMTME 

and COMNME. 
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TABLE B4.- SUBROUTINE COMMON REQUIREMENTS 


ROUTINE 

*DECK 

COMTME 

COMVAR 

ACSIM 

ACMOD 

X 

X 

AERO 

START 


X 

ANGL 

START 


X 

ARIN 

DEBUG 


X 

ARNM 

DEBUG 

X 

X 

ATRK 

START 



AVAR 

MISCL 



CINS 

START 


X 

CONV 

MISCL 


X 

COOR 

START 


X 

EVALF 

EVLMD 


X 

EVALH 

EVLMD 


X 

EVLFW 

EVLJC 


X 

EVLFX 

EVLJC 


X 

EVLHX 

EVLJC 


X 

FILT 

MISCL 



FOIL 

START 


X 

INIT 

SETUP 

X 

X 

JACSX 

SETUP 



MINC 

OPTIM 



MODL 

ACMOD 

X 

X 

NOYZ 

ACMOD 


PGPLO 

PLOTS 



PLOT 

PLOTS 

X 

X 

PRNT 

PRNTS 


X 

RADR 

START 


X 

READ 

SETUP 

X 

X 

SECOND 

MISCL 



SENS 

OPTIM 

X 

X 

SETP 

OPTIM 

X 

X 

SOLU 

PRNTS 

X 

X 

STRT 

START 

X 

X 

TRAJ 

OPTIM 

X 

X 

TRANS 

EVLMD 


X 

VELS 

START 


X 

VRNC 

PRNTS 

X 

X 

WIND 

START 


X 

XYPLO 

PLOTS 




COMNME 


X 

X 


X 

X 

X 

X 

X 

X 

X 

X 


X 

X 


Jacobian Calculations 

Calculations involving the Jacobian matrices defined in chapter 2 are a large 
part of the computational burden for the SMACK algorithm. The Jacob ians are sparse 
and special techniques are employed in SMACK to make their use as efficient as 
possible. If a user should wish to make a change in the state or output model he 
will need to know how to modify the corresponding Jacobian arrays. This section 


91 


gives the details necessary to understand the Jacobian structure. All calculations 
for the backward filter and forward smoother are carried out in subroutine SENS, the 
Jacobian arrays FX, FW, and HX are evaluated in EVLFX, EVLFW, and EVLHX, respec- 
tively, which are called by SENS. 

The Jacobians and their companion arrays are found in COMMON block /JAC/ 

(table B2 ) . A description of Jacobian sets (FX, KFX, LFX, MFX), (FW, KFW, LFW, 

MFW) and (HX, KHX , LHX , MHX) will be illustrated by considering the first set 
only: the other two sets are similarly structured. The key for the row-column 

location for any element of FX is given in table B5. The corresponding row index is 
stored in LFX, the column index in MFX. For example, the fifth element o , w ic 
corresponds to the fifth partial-derivative expression to be evaluated in EVLFX is, 
in fact, the partial derivative of Y1DT with respect to Y2. This is seen to be the 
fourteenth row and fifth column of the Jacobian matrix FX. Hence, 


LFX(5)=14; MFX(5)=5 


TABLE B5.- KEY FOR LOCATION OF FX ENTRIES 
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These entries are verified by inspection of the DATA statements for arrays LFX 
and MFX found in subroutine EVLFX and repeated here: 


DATA ( LFX (I), 1=1, 72 )/10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 1, 2, 

3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15,13, 
!4, 15, 13, 14, 15, 13, 14, 15, 13, 14, 15, 13, 14, 15, 

* 13,14,15,13,14,15,13,14,15,13,14,16,16,17, 

18, T6, 17, 18, 16, 16, 17, 18, 16, 17, 18, 16, 17, 18, 

* 16 , 18 / 

DATA ( MFX ( I ) , I = 1 , 72)/ 1, 2, 3, 4, 5, 6,10,11,12,13,14,15,22,23, 

24,25,26,27,28,29,30,31,32,33,34,35,36,52, 
52,52,53,53,53,54,54,54,82,82,82,83,83,83, 
84,84,84,16,16,16,17,17,17,18,18,55,56,56, 

56,57,57,57,85,86,86,86,87,87,87, 16,16,16, 
17,17/ 


Tables B6 and B7 give the corresponding structural details for the FW and HX 
Jacobians. The DATA statements from subroutine EVLFW are 


DATA(LFW(I), 1=1,15)/ 1, 
DATA(MFW(I),I=1,15)/ 1, 


2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15/ 
2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15/ 


and the DATA statements from subroutine EVLHX are 


DATA(LHXd), 1.-1,136)/ 1, 2, 3, 4, 5, 6, 7, 7, 7, 8, 8, 9, 9, 9, 

10.10.11.11.12.13.13.13.14.14.14.15.15.15, 

13.13.13.14.14.14.15.15.15.13.13.13.14.14, 

14.15.15.15.13.13.13.14.14.14.15.15.15, 0, 

16, 16, 17,17, 17, 18, 18,18,16,16,16,17,17, 17, 

18.18.18.16.16.16.17.17.17.18.18.18.16.16, 

16.17.17.17.18.18.18.19.19.20.20.21.21.19, 
20,20,21,21,22,22,23,23,24,24,22,22,23,23, 
2 3,24,24,24,22,23,23,24,24,25,25,25,26,26, 
27,27,27,28,28,28,29,29,30,31/ 

DATA(MHX( I ) , 1= 1 , 136 ) / 16 , 17 , 18, 1 9 ,20, 2 1 , 1 9 ,20,21 , 19,20, 19 20, 2 1 

7, 8, 7, 8, 9, 7, 8, 9, 7, 8, 9, 7, 8, 9) 

10.11.12.10.11.12.10.11.12.13.14.15.13.14, 
15,13,14,15,16,17,18,16,17,18,16,17,18, 0, 
2 ’ 3 > !» 2 > 3, 1, 2, 3,10,11,12,10,11,12, 
10,11,12, 4, 5, 6, 4, 5, 6, 4, 5, 6,16,17, 

18.16.17.18.16.17.18.10.12.11.12.11.12.17, 

16.17.16.17, 1, 3, 2, 3, 2, 3,11,12,10,11, 

12.10.11.12.17.16.17.16.17.19.20.21.19.20, 
19,20,21,16,17,18,13,14,13,14/ 


Note that location 56 is not used, 
are provided in subroutine INIT. 


Data for the VB and SF locations (137-196) of HX 
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TABLE B6 . - KEY FOR LOCATION OF FW ENTRIES 



ODD 
L M N 

D D D 
X Y H 

G 6 G 
X Y H 

P T P 
H H S 
2 2 2 

X Y H 
2 2 2 


PH2T 
TH2D 
PS 20 

1 

2 

3 





1 

2 

X2DT 

Y2DT 

H2DT 


4 

5 

6 




4 

5 

6. 

WXDT 

WYDT 

WHDT 



7 

8 

9 



7 

8 
9 

PH1D 

TH1D 

PS1D 




10 

11 

1 L 


10 

11 

12 

X1DT 

Y1DT 

H1DT 





13 

14 

15 

13 

14 

15 

PHID 

THTD 

PSIO 

XDOT 

YDOT 

HOOT 








1 2 3 

4 5 6 

7 8 9 

10 11 12 

13 14 15 



It should be emphasized that the array structures shown in tables B5-B7 are 
master lists: any particular problem will use subsets of the partial-derivative 

expressions in subroutines EVLFX, EVLFW, and EVLHX. The selection of those to be 
included is done in subroutine READ, following the coding list analysis. The index 
arrays KFX , KFW, and KHX are used for this purpose: the location of ^ Partia 
derivative to be evaluated is set equal to the array location of FX, FW or HX that 
it will fill or else the location is set to zero. For example, the fifth partial 
derivative expression in EVLFX is set up as 


K=KFX(5) 

IF(K.GT.O) FX(K) = 1 

In this way the Jacobian subset in use completely fills the first NFX locations of 
FX Corresponding elements of LFX and MFX are modified accordingly. Examples of 
the Jacobian setup are shown with the output for the test problem given in 

appendix C. 
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TABLE B7 . - KEY FOR LOCATION OF HX ENTRIES 



P T P 
H H S 
2 2 2 

X Y H 
2 2 2 

W W V 
X Y h 

1 P T F 
H H F 
1 1 1 

> X Y 
i 1 1 

H P T F 
H H S 
I T I 

X Y H 

! 

V V 

B to B 
1 30 

s s 

F to F 
1 30 

■ i 

PHI 

THT 

PSI 






1 

2 

3 


137 

138 

139 

167 

168 
169 

1 

2 

3 

X 

Y 

H 







4 

5 

6 

140 

141 

142 

170 

171 

172 

4 

5 

6 

RNg 

BRG 

ELV 







7 8 9 
10 11 
12 13 14 

143 

144 

145 

173 

174 

175 

7 

8 
9 

WXY 

WHO 

VWD 



i5 16 
17 18 

19 





146 

147 

148 

176 

177 

178 

10" 

11 

12 

VT 

AV 

BV 



20 21 22 
23 24 25 
26 27 28 

29 30 31 
32 33 34 
35 36 37 

ib 39 40 
41 42 43 
44 45 46 

47 48 49 
50 51 52 
53 54 55 


149 

150 

151 

179 

180 
181 

13 

14 

15 

Ax 

AY 
A Z 

57 58 
[ 59 60 61 
62 63 64 

74 75 76 
77 78 79 
80 81 82 


[55 66 67 
68 69 70 
71 72 73 


83 84 85 
86 87 88 
89 90 91 


"152 

153 

154 

~ rer 

| 183 

184 

nr* 

17 

18 

P 

Q 

R 




~52 93 

94 95 
96 97 


6H 

99100 

101102 


155 

156 

157 

TBS 

186 

187 

T3H 

20 

21 

AL 

AM 

AN 

105 104 

105106 
107108 



T05TT8 

11112113 

14115116 


— TT7 — 

118119 

120121 


133 

159 

160 

T58 

189 

190 

55 r 

23 

24 

RNA 

BRA 

ELA 







72123124 

25126 

27128129 

“161 

162 

163 

— 

192 

193 

55" 

26 

27 

HOG 

VGR 

TRK 




1 

1 

33134 

35136 

30131132 


'164 

165 

166 

"184 

195 

196 

29 

30 


1 2 3 

4 5 6 

7 8 9 

10 11 12 

13 14 15 

16 17 18 

19 20 21 ; 

57 - 66 i 

67 - 96 



Output Modification 

The mathematical representation for an aircraft utilized in SMACK may not be 
sufficient for every application. Since the measurement model is the most likely 
area to require modification, an example is considered here that illustrates the 
addition of a new measurement to the set discussed in chapter 3 and listed in 
table 5.1. The new measurement is to D e the heading from a directional gyro. This 
instrument is a two-degree-of-freedom gyro, and has its outer gimbal axis parallel 
to the body vertical axis, with its spin axis aligned with the magnetic North 
Pole. The heading, as measured by the outer gimbal angle, is a function of the 
aircraft attitude (roll, pitch, and yaw). The output relation is summarized in 
table B8, which also includes expressions for partial derivatives of heading with 
respect to the Euler angles. 


TABLE B8 . - TWO-DEGREE-OF-FREEDOM DIRECTIONAL GYRO 


Output Relation 

!| . = tan -1 [(cos* sin* m - sin* sine cos* m /cos0 cos* m ] 

O 


where 


* m = * + A* 

and (*, 0, *) are the Euler angles, * m is the magnetic heading, and A* is the 
magnetic variation. The user must specify the magnetic variation by setting 
parameter VARM in subroutine DATA (see appendix E). 


0* g /3*) 
( 3* g /ae) 

(3*g/3*) 


Partial Derivatives 

-cos0 cos* m (cos* sin© cos* m + sin* sin* m )/D 
cos* m (cos* sinO sin* m - sin* cos* m )/D 
cos* cos0/D 


where 

P 2 

D = (cos* sin* m - sin* sin0 cos* m ) + (cos0 cos* m ) 


The changes made to the program were relatively simple: only COMMON group 

COMNME and subroutines READ, INIT, STRT, EVALH, and EVLHX required modification. 

The name HDG was chosen to designate the new variable; it now occupies the formerly 
unused 28th location of the output vector (element 52 of the variable arrays). This 
choice required that HDG be equivalenced to ED(52) in COMNME (see table B3), and 
that the DATA statement for the NME array in subroutine READ include 

NME(52)=3HHDG 

and the conversion factor (0.01745) and unit name (DEG) be assigned in the DATA 
statements for the ISCE and IUNE arrays (in BLOCK DATA, following INIT) by setting 

ISCE(52)=3 ; IUNE(52)=7 

Note that these integers point to locations of arrays SKL and UNM (in /DEF/) 
defined in BLOCK DATA as 
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C/2- 


SKL(3)=0. 01745 ; UNT( 7 ) =4HDEG 



Now, to accommodate the requirements for a starting solution should there be no 
direct measurement of yaw angle (PSI), there has been inserted in subroutine READ 
the following code: 

C SET UP OUTPUT ARRAY FOR HEADING ESTIMATE 

IF(KH(28).EQ. 1) KH(3)=1 

and in subroutine STRT the code: 

I F( KZ ( 3 ) .NE. 1 .AND.KZ(28) .NE.O) THEN 

C PUT HDG ESTIMATES INTO PSI CHANNELS FOR STARTING SOLUTION 

DO 80 N= 1 , NPTS 
H(3,N)=H(28,N) 

D(3,N)=D(28,N) 

80 CONTINUE 
KZ(3)=-1 
END IF 

Next, the following statements defining the output relation defined in table B8 were 
added to subroutine EVALH: 

IF(KH(28) .GT.O) THEN 

C PERFORM TWO-DOF GYRO CALCULATIONS 

PSM=PSI+VARM 

SNPSM=SIN(PSM) 

CSPSM=C0S ( PSM ) 

A 1 =CSPHI*SNPSM-SNPHI*SNTHT*CSPSM 

A2=CSTHT*CSPSM 

HDG=ATAN2( A 1 , A2 ) 

END IF 

Finally, the three new elements of the output Jacobian (HX, KHX , LHX , MHX) were 
assigned to locations 130-132 (see table B7) and provision for evaluation of the 
partial derivatives from table B8 is accomplished by inserting into subroutine EVLHX 
the following statements: 
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IF(KH(28).GT.O) THEN 


PSM=PSI+VARM 
SNPSM=SIN( PSM) 

CSPSM=COS( PSM) 

A1 =CSPHI*SNPSM-SPHST*CSPSM 
A2=CSTHT*CSPSM 
DSQ=A1 **2+A2**2 
CHDG=SF(28)/DSQ 


C PARTIALS OF HDG WRT PHI, THT, PS I 


K=KHX( 130) 
IF(K.GT.O) 
K=KHX( 131 ) 
IF(K.GT.O) 
K=KHX( 1 32 ) 
IF(K.GT.O) 


HX( K ) =CHDG* ( - A2 ) * ( CPHST*CSPSI+SNPHI*SNPSI ) 
HX( K) =CHDG*CSPSI*( CPHST*SNPSI-SNPHI*CSPSI ) 
HX( K) =CHDG*CPHCT 


END IF 

Use of the heading measurement in an accident situation is illustrated by the 
second problem of appendix D. There, records of airspeed, altitude, heading, and 
vertical acceleration recovered from an aircraft flight data recorder are combined 
with radar tracking data to estimate forces and motions along the accident 
trajectory. 


Temporary Storage 

When the SMACK algorithm of table 2.1(b) estimates forcing-function time his- 
tories, the arrays d(i) and L(i), calculated during the backward-filter pass, must 
be stored for use during the forward-smoother pass. These steps of the algorithm 
are performed in subroutine SENS, where d(i), an NW*1 vector, and L(i), an NW*NXP 
matrix are combined in a vector QFWTBP of length NW*(NXP+1). Here i refers to the 
time point, NW is the number of forcing functions (NW < 9), and NXP is the number of 
parameters (NXP < 42). Hence, a maximum of 387 words must be stored at each of a 
possible 6000 time points, requiring a total of 2,322,000 words of temporary 
storage. 

The method of acquiring large amounts of high-speed scratch storage is usually 
machine dependent. Table B9 illustrates three ways of realizing and accessing the 
temporary storage needed in SENS. The first method shown is probably limited to 
systems that have virtual memory. It has been used in the VAX versions of SMACK. 

The second and third methods are suitable for use with either VAX or Cray ver- 
sions. The storage required is provided on the Cray system by high-speed, solid- 
state memory. The third method shown in table B9 is the most efficient, and is now 
used in both Cray and VAX versions of SMACK. The statements that open the scratch 
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TABLE B9 . - REALIZATION OF TEMPORARY STORAGE 


Method 1 

DIMENSION SCRATCH( 387, 6000) 
DIMENSION QFWTBP(387) 


DO 100 L = 1 , NWRD 
SCRATCH(L,N)=QFWTBP(L) 
100 CONTINUE 


DO 200 L=1,NWRD 
QFWTBP ( L ) =SCRATCH ( L , N ) 
200 CONTINUE 


Method 2 

DIMENSION QFWTBP(387) 


WRITE(3) (QFWTBP(L) ,L=1 ,NWRD) 


BACKSPACE(3 ) 

READ( 3) (QFWTBP(L) ,L=1 ,NWRD) 
BACKSPACE( 3) 


Method 3 

DIMENSION QFWTBP (387) 

WRITE( 3 , REC=N) (QFWTBP(L) ,L=1 ,NWRD) 
RE AD ( 3 , REC=N ) ( QFWTBP (L ) , L= 1 , NWRD ) 
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file are found in the main program: 


C OPEN SCRATCH FILE USED IN SR SENS (BF-FS) 

NBYT=8 !CRAY/VAX DP SPECIFIC 
LREC=NWT*(NXP+1 )*NBYT 

0PEN(UNIT=3,STATUS= 'SCRATCH' , ACCESS= ' DIRECT' , RECL=LREC ) 

where LREC is the record length in bytes, and NBYT is the number of bytes per 
word. Note that the number of bytes per word is machine dependent, i.e., for the 
VAX single-precision version, NBYT=4 , while for the Cray and VAX double-precision 

versions, NBYT=8. 
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APPENDIX C 


FLIGHT-TEST PROBLEM 


A test problem is presented here to further illustrate the important features 
of a solution by SMACK, and to provide a reference output listing for checking a new 
installation of the program. The internal aircraft-trajectory simulation (a rising, 
coordinated 180 turn in winds) written into subroutine MODL is utilized to provide 
the data record. The problem solved is typical of a data-consistency analysis 
ollowing a flight-test experiment (as described in chapter 6). Radar tracking of 
aircraft position is included in the measurements, thus making it possible to 
estimate winds along the flightpath. The description of the solution and its output 
will follow the steps outlined in the block diagram shown in figure A2. The 
examples of SMACK output given here and in appendix D were run on the Cray system. 

A coding list for the test problem is shown in figure Cl. In addition to the 
radar tracking data (RNG, BRG, ELV), there are measurement records available from 

?PHT ar ?HT nS p^n ent ! ?D Cl n din ? accelerometers ( AX > AY, AZ), attitude and rate gyros 
(PH 1 , THT, PSI) and (P, Q, R), altimeter (H), and aerodynamic angle vanes (AV 

BV). The record of dynamic pressure is assumed to have been converted to true 
airspee (VT). For this problem, all the data records are to be fit (1=1 J=1 on 
the VAR statement), except for the rate-gyro measurements, which are to be used as 
orcmg functions (1 = 1, J=0 on the VAR statement). A number of bias errors and 

(GX ry aC ^r b nv ( dete ™ ined ’ along with time histories of forcing functions 

ihp’pnH 1 GH . a " h (DX ’ ’ ' M ° te that the S6COnd gr ° up need nofc be mentioned in 

e coding list because accelerometer measurements (AX, AY, AZ) are included in the 

data set. Note also that the low-pass filter cutoff frequency for roll angle (PHI) 

^he sa e :n 1 SPeCif i ed * * °:1 5 ^ ^ "° minal Value is °- 10 (as calculated from 

the sampling interval specified on the REC statement). The nominal value will be 

used to fUtor all other records. Choice of a value different than nominal is made 

MKS anttn Tr \ White " residual f°r construction of a starting solution. The 
MKS and END statements indicate that the solution option flags will be set as 

ISTRT=0, IS0LU=0, IDATA= 1 , IPL0T=0, IDBUG=1 

A coding list summary is printed by subroutine READ, as shown in fig- 
ure C1(b) When READ has finished setting up all the integer arrays, control is 
transferred to subroutine ARIN, which prints them. The first group of these the 
state-variable, forcing-function, and output-variable integer arrays, are shown in 
figure C2 An interpretation of the output-variable arrays will be given here 
because these are directly determined by coding-list entries. The first column of 
figure C2(c) refers to the order of output variables listed in table 5 1 The 
second and third columns are the arrays that indicate which variables are to be 

me^hp’tot ? ent K y in r the KH a " ay iS ° ne if J=1 in the corresponding VAR state- 
ent The total number of output variables to be estimated is seen to be NH=18 

The IH array contains a (compressed) list of addresses for all variables to be 
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(b) 

SUMMARY 

OF 

CODING 

LIST: 



NF - 15 (DIFFERENTIAL EQUATIONS IN STATE MODEL) 
NX - 15 (ESTIMATED VARIABLES IN STATE MODEL) 

NW - 6 (ESTIMATED FORCING FUNCTIONS) 

NH - 18 (ESTIMATED VARIABLES IN OUTPUT MODEL) 

NZ - 16 (AVAILABLE MEASUREMENTS) 

NV - 13 (OUTPUT ESTIMATES FIT TO MEASUREMENTS) 
NZH- 3 (OUTPUTS MEASURED. NOT ESTIMATED) 

NHZ- 21 (NH + NZH) 

NZX- 3 (STATES MEASURED. NOT ESTIMATED) 

NXZ- 18 (NX + NZX) 

NXO- 15 (STATE-VARIABLE INITIAL CONDITIONS) 
NWB- 0 (FORCING-FUNCTION MEAN VALUES) 

NVB- 6 (INSTRUMENT BIAS ERRORS) 

NSF- 3 (INSTRUMENT SCALE FACTORS) 

NT - 24 (NXO + NWB + NVB + NSF) 


Figure Cl.- Flight-test problem, (a) Coding list, (b) summary. 

estimated. The next two columns, arrays KZ and IZ, similarly identify those vari- 
ables specified as measured in the coding list (an entry in array KZ is one if 
1 = 1 on P the VAR statement). The sixth column, array IV, identifies those vanab es 
that are both measured and estimated. Hence, for this problem, NV=13 estimates 
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(a) STATE VARIABLE ARRAYS: 
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Figure C2.- Contents of integer arrays for flight-test problem, (a) State 

variables, (b) forcing functions. 

will be "fit" to the measurement records during the solution. The next column 
(array IZH) lists those variables measured, but not to be estimated (NZH=3) and the 
following column (array IHZ) lists all output variables mentioned in the coding list 


103 


(c) OUTPUT VARIABLE ARRAYS: 

I KH IH KZ IZ IV IZH IHZ KVB IVB KSF ISF 

1 i i l 1 119 10 16 0 19 

212122 20 20 17 0 20 

3 1 3 1 3 3 21 30 18 0 21 

414066040 19 00 

515077050 20 00 

616188060 21 00 

717199070000 
8 1 8 1 13 13 080000 

9191 14 14 090000 

10 1 10 0 15 15 0 10 0 0 0 0 

11 1 11 0 16 16 0 11 0 0 0 0 

12 1 12 0 17 17 0 12 0 0 0 0 

13 1 13 1 18 18 0 13 0 0 0 0 

14 1 14 1 19 0 0 14 0 0 0 0 

15 1 15 l 20 0 0 15 0 0 0 0 

16 1 16 1 21 0 0 16 1 0 0 0 

17 1 17 1 000 17 1 000 

18 1 18 1000 18 1000 

19 001000 19 1010 

20 001000 20 1010 

21 001000 21 1010 

22 00000000000 

23 00000000000 

24 00000000000 

25 00000000000 

26 00000000000 

27 00000000000 

28 00000000000 

29 00000000000 

30 00000000000 

Figure C2.- Concluded, (c) Output variables. 

(NHZ=21). Information about the bias-error and scale-factor parameters to be iden- 
tified (K=1 and/or L=1 on VAR statement) is stored in the arrays shown in the last 
four columns. Note that NVB=6 and NSF=3 for this problem. 

The next section of output displayed by subroutine ARIN is shown in fig- 
ure C3. There the contents of the arrays defining the structure of the Jacobians 
FX, FW, and HX are listed. The method for Jacobian representation was discussed 
earlier in appendix B: the Jacobian matrices are usually quite sparse, and it is 
convenient to store the nonzero elements of each in a linear array. Here the method 
is illustrated by considering the representation of HX, the output Jacobian. Hence, 
although HX in figure C3(c) is referred to as an NV*NXP matrix, it is actually 
represented in the program by an NHX-element vector . Addresses of the solution- 
active elements of HX are stored in the linear array KHX, as shown in figure C3(c). 
where consecutive locations are listed from left to right. Corresponding row and 
column addresses for each of the NHX active elements of HX are stored in the linear 
arrays LHX and MHX, also listed from left to right. This information is summarized 
in the NV*NXP array that represents the actual row-column locations of consecutive 
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(a) JACOBIAN FX(NF*NT): NF- 15 , NT- 24 , NFX- 25 , NFXT- 72 
ARRAY KFX (SPECIFIES ACTIVE ELEMENTS) 

000123000456000000000 

000000000000000000000 

0 0 0 0 0 0 0 0 0 0 0 7 8 9 10 11 12 13 14 15 16 

17 18 19 20 21 22 23 24 25 

ARRAY LFX (ROW INDEX) 

7 8 9 13 14 15 10 10 11 12 10 11 12 10 10 11 12 10 11 12 10 

11 12 10 12 

ARRAY MFX (COL INDEX) 

1 2 3 7 8 9 19 20 20 20 21 21 21 22 23 23 23 24 24 24 10 

10 10 11 11 


X2DT 

Y2DT 

H2DT 

WXDT 

WYDT 

WHDT 

X1DT 

Y1DT 

H1DT 

PHID 

THTD 

PSID 

XDOT 

YDOT 

HDOT 


ARRAY FX(L, M) - FX(KFX(K)), KFX(K) > 0 


X 
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00000006000 


XYHVVVVVVSSS 
BBBBBBFFF 
16 17 18 19 20 21 19 20 21 

000000000000 
000000000000 
000000000000 
000000000000 
000000000000 
000000000000 
000000000000 
000000000000 
000000000000 
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00000009 12 0 16 19 

0 0 0 0 0 0 0 10 13 0 17 20 

000000000000 
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000000000000 


Figure C3. Contents of Jacobian arrays for flight-test problem. 

(a) State model FX. 


entries of HX. Note that the NXP parameters defining the columns comprise 15 state 
initial conditions, 6 bias errors, and 3 scale factors. The NV rows of the array 
correspond to the output variables that are to be fit to corresponding measurement 
records. Note also that because the rate-gyro measurements are used as forcing 
functions in this problem, the partial derivatives of P, Q, and R with respect to 

their respective bias errors and scale factors are associated with the state model 
not the output model . ' 
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(b) JACOBIAN FW(NF*NW): NF= 15 , NW= 6 , NFW- 6 , NFWT- 15 
ARRAY KFW (SPECIFIES ACTIVE ELEMENTS) 
000123456000000 
ARRAY LFW (ROW INDEX) 

1 2 3 4 5 6 

ARRAY MFW (COL INDEX) 

1 2 3 4 5 6 

ARRAY FW(L,M) - FW(KFW(K)), KFW(K) > 0 

D D D G G G 

X Y H X Y H 


X2DT 1 0 0 0 0 0 

Y2DT 0 2 0 0 0 0 

H2DT 0 0 3 0 0 0 

WXDT 0 0 0 4 0 0 

WYDT 0 0 0 0 5 0 

WHDT 0 0 0 0 0 6 

X1DT 0 0 0 0 0 0 

Y1DT 0 0 0 0 0 0 

H1DT 0 0 0 0 0 0 

PHID 0 0 0 0 0 0 

THTD 0 0 0 0 0 0 

PSID 0 0 0 0 0 0 

XDOT 0 0 0 0 0 0 

YDOT 0 0 0 0 0 0 

HDOT 0 0 0 0 0 0 

Figure C3.- Continued, (b) State model FW. 

The output summarizing the trajectory for the simulated aircraft maneuver 
provided by subroutine MODL for the measurement record is shown in figure C4. The 
initial-condition set for the active elements of the state model is shown first, 
followed by bias-error and scale-factor values. Next are listed the statistics of 
random noise added to the output variables to simulate measurement errors. Each 
value in the STDEV column (specified by VAL1 on the VAR statement) is used to form 
the corresponding diagonal element of the measurement-error covariance matrix. Note 
that each measurement is biased by the value shown in the column headed MEAN. The 
mean and rms values of each unmeasured forcing function are also computed and 
listed. The rms value for each of these variables is used to form the corresponding 
diagonal element of the unknown forcing-function covariance matrix. The final data- 
record information gives the number of points per record, the length of the record, 


106 



(c) JACOBIAN HX(NV*NT) : NV= 13 , NT= 24 , NHX- 60 , NHXT-192 


ARRAY KHX (SPECIFIES ACTIVE ELEMENTS) 


1230045678 
15 16 17 18 19 20 21 0 0 0 
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0000000000 
51 52 53 54 55 56 57 0 0 0 
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Figure C3.- Concluded, (c) Measurement model HX. 


and the nominal bandwidth for the low-pass filter used to 
in the STRT subroutine. 


process the measurements 
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DATA FROM AN INTERNAL SIMULATION 

THE TRAJECTORY IS AN A/C MANEUVER CONSISTING 
OF A RISING COORDINATED 180-DEG TURN IN WIND 


INITIAL CONDITIONS: 


VARIABLE 


VALUE 


STDEV 


X2 , 

M/S2 

0 . 5410E-02 


0 . OOOOE+OO 

Y2 , 

M/S2 - 

0. 1534E-02 


O.OOOOE+OO 

H2 , 

M/S2 

0 . 3206E-03 


O.OOOOE+OO 

WX , 

M/S 

0 . 2081E+01 

+ , “ 

O.OOOOE+OO 

WY , 

M/S 

-0 . 4546E+01 

+ , “ 

O.OOOOE+OO 

WH , 

M/S 

0. 1000E+01 


O.OOOOE+OO 

XI , 

M/S 

O.OOOOE+OO 

+ i “ 

O.OOOOE+OO 

Y1 , 

M/S 

0 . 8500E+02 

+ * “ 

O.OOOOE+OO 

HI , 

M/S 

O.OOOOE+OO 

+ •“ 

O.OOOOE+OO 

PHI, 

DEG 

O.OOOOE+OO 

+ ,~ 

O.OOOOE+OO 

THT , 

DEG 

0. 1422E+01 

+ » " 

O.OOOOE+OO 

PSI, 

DEG 

0 . 9135E+02 


O.OOOOE+OO 

X . 

NM 

-0 . 9112E+00 


O.OOOOE+OO 

Y , 

NM 

-0 . 5400E+00 

+ •- 

0 . OOOOE+OO 

H , 

M 

0. 1000E+04 

+ .~ 

0 . OOOOE+OO 

PHI, 

D/S 

-0 . 9415E-01 

+ > ” 

O.OOOOE+OO 

TH1 , 

D/S 

0 . 2079E-03 


0 . OOOOE+OO 

PSI, 

D/S 

-0. 1163E-01 


O.OOOOE+OO 

BIAS 

ERRORS 

: 



VARIABLE 

VALUE 


STDEV 

AX , 

G'S 

O.OOOOE+OO 


O.OOOOE+OO 

AY , 

G'S 

O.OOOOE+OO 

+ •- 

O.OOOOE+OO 

AZ , 

G'S 

O.OOOOE+OO 


O.OOOOE+OO 

P , 

D/S 

O.OOOOE+OO 

+ * “ 

O.OOOOE+OO 

Q , 

D/S 

O.OOOOE+OO 


O.OOOOE+OO 

R , 

D/S 

O.OOOOE+OO 

+ i ~ 

O.OOOOE+OO 


SCALE FACTORS: 
VARIABLE VALUE 


STDEV 


P , D/S 
Q , D/S 
Q , D/S 
R , D/S 


0 . 1000E+01 +,- 0 . OOOOE+OO 
0 . 1000E+01 +,- 0. 0000E+00 
0 . 1000E+01 +,- O.OOOOE+OO 
0 . 1000E+01 +,- O.OOOOE+OO 


Figure C4.- Output summarizing simulated flight-test maneuver. 


After acquiring the data record, SMACK passes program control to subroutine 
STRT which determines a set of initial conditions and forcing-function time his- 
tories necessary to generate a starting solution for the SMACK algorithm. The block 
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MEASUREMENT NOISE: 


VARIABLE 

MEAN 


STDEV 

PHI, 

DEG 

-0 . 3289E-02 

+ . - 

0 . 5000E-01 

THT, 

DEG 

-0 . 2518E-02 

+ , - 

0. 5000E-01 

PSI, 

DEG 

-0. 1502E-02 


0 . 5000E-01 

H , 

M 

0. 1099E-01 


0. 1000E+00 

RNG , 

NM 

-0. 1634E-03 

+ » - 

0 . 5000E-02 

BRG, 

DEG 

-0. 1458E-02 


0 . 5000E-01 

ELV, 

DEG 

0 . 6007E-03 


0 . 5000E-01 

VT , 

KT 

0 . 4505E-02 

+ , ~ 

0. 1000E+00 

AV , 

DEG 

-0. 1008E-01 


0. 5000E-01 

BV , 

DEG 

0 . 5341E-04 


0 . 5000E-01 

AX , 

G'S 

-0. 1700E-03 


0. 1000E-01 

AY , 

G'S 

0 . 6972E-03 


0. 1000E-01 

AZ , 

G'S 

0 . 4576E-03 


0. 1000E-01 

P , 

D/S 

0. 1719E-03 


0. 1000E-02 

<? . 

D/S 

0 . 8182E-04 

+ . - 

0. 1000E-02 

® > 

D/S 

-0. 1006E-03 

+ . - 

0. 1000E-02 


FORCING FUNCTIONS (UNMEASURED): 


VARIABLE MEAN 


STDEV 


DZ , M/S3 
DY , M/ S3 
DH , M/S3 
GX , M/S2 
GY , M/S2 
GH , M/S2 


-0. 1202E-03 
-0 . 5255E-11 
-0 . 7124E-05 
-0. 1579E-15 
-0 . 2368E-15 
0.7895E-16 


+,- 0 . 1643E+00 
+,- 0 . 9321E-01 
+,- 0 . 9736E-02 
+,- 0 . 8128E-01 
+.- 0. 9604E-01 
+,- 0.4936E-01 


DATA RECORD INFORMATION: 


90 ( LENGTH OF RECORD, POINTS) 
®?£L- 0 • 9000E+02 (LENGTH OF RECORD, SECONDS) 
FREQ- 0 . 1000E+00 (LP FILTER BANDWIDTH, HZ) 


Figure C4.- Concluded. 


diagram describing this procedure was shown in figure A3. Generally, each measure- 
ment record is filtered with the cutoff frequency set to make the residual appear 
white. The residuals are used to determine the diagonal elements of R the 
measurement-error weighting matrix. For a solution using simulated data , ’however 
the weights used are the actual measurement-noise covariance values The 
filtered records are then used to estimate a set of initial conditions for the 
trajectory, and to construct a set of forcing-function time histories. The variance 
of each is used as the appropriate diagonal element of weighting matrix Q, unless 
the data has been simulated, in which case the "true" value is used. Results of 
starting solution calculations are shown in figure C5. 
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STARTING SOLUTION TABULATION 
RESPONSE ERRORS AND WEIGHTS : 


VARIABLE MEAN STDEV 


PHI, 

DEG 

0 . 6384E-13 

+ , ” 

0 . 4408E-01 

THT , 

DEG 

0. 1598E-13 

+ » - 

0 . 4442E-01 

PSI, 

DEG 

-0. 1949E-12 

+ » “ 

0 . 5192E-01 

H 

M 

-0 . 4689E-11 

+ * " 

0 . 8529E-01 

RNG, 

NM 

0. 5893E-14 

+ , " 

0 . 4173E-02 

BRG, 

DEG 

0 . 2172E-11 

+ > ~ 

0 . 5154E-01 

ELV, 

DEG 

-0 . 1697E-13 

+ , “ 

0 . 4268E-01 

VT , 

KT 

0 . 9823E-13 


0. 1070E+00 

AV , 

DEG 

-0 . 4468E-13 


0 . 4081E-01 

BV , 

DEG 

-0 . 7069E-16 


0 . 4509E-01 

AX , 

G'S 

0. 1540E-15 


0 . 8203E-02 

AY , 

G'S 

0. OOOOE+OO 


0 . 8688E-02 

AZ , 

G'S 

0 . 6440E-14 

+ * ” 

0 . 8911E-02 


FORCING FUNCTIONS AND WEIGHTS: 


VARIABLE MEAN 


STDEV 


DX , M/ S3 

0 . 9814E-03 

+ > ~ 

DY , M/S3 

-0 . 3048E-03 


DH , M/ S3 

-0 . 3726E-03 


GX , M/S2 

-0. 5356E-01 

+ , “ 

GY , M/S2 

-0. 1677E-01 


GH , M/S2 

-0. 1078E-02 

+ , ” 

PHI, D/S 

-0 . 8461E-02 

+ >" 

TH1 , D/S 

0 . 2532E-01 

+ > “ 

PSI, D/S 

- 0 . 2033E+01 

+ » “ 


0. 1639E+00 
0 . 9270E-01 
0. 1599E-01 
0 . 8816E-01 
0 . 1050E+00 
0. 5965E-01 
0 . 7375E+00 
0. 1357E+00 
0. 9920E+00 


Figure C5.- Tabulation of starting solution 


(R* *0 . 5) 

0 . 5000E-01 
0 . 5000E-01 
0. 5000E-01 
0 . 1000E+00 
0 . 5000E-02 
0 . 5000E-01 
0 . 5000E-01 
0. 1000E+00 
0 . 5000E-01 
0 . 5000E-01 
0. 1000E-01 
0 . 1000E-01 
0. 1000E-01 


(Q**0.5) 

0. 1643E+00 
0 . 9321E-01 
0.9736E-02 
0 . 8128E-01 
0 . 9604E-01 
0 . 4936E-01 
0 . 0000E+00 
0 . OOOOE+OO 
0. 0000E+00 


FILTER 

0 . 1500E+00 
0 . 1000E+00 
0 . 1000E+00 
0 . 1000E+00 
0. 1000E+00 
0 . 1000E+00 
0 . 1000E+00 
0. 1000E+00 
0. 1000E+00 
0. 1000E+00 
0. 1000E+00 
0. 1000E+00 
0. 1000E+00 


for flight-test problem. 


Following the starting solution, the program calls subroutine ARNM which loads 
initial state-variable and forcing-function values into variable arrays XD, WD 
displays them, and calls EVALF, which evaluates the state derivative FD, and t en 
calls EVLFX and EVLFW , which evaluate the state-model Jacob lans FX and FW, and 
displays the contents of these arrays. The purpose of this output is for debugging 
and/or general information only. The contents of arrays XD, WD, FD FX, and FW are 
shown in figure C6(a), where the element order is from left to right across the 
page. The variables are referenced in tables 5.2, 5.3, and B4. Units are the 
internal (mks-radian) units used by SMACK. The first integration step is t en 
performed, and the results used to evaluate the output model and output Jacobian y 
calling EVALH and EVLHX . The output of arrays HD, SD, and HX is shown in lg 
ure C6(b). These variables are referenced in tables 5-1, B3, and B . 

The problem is now ready for solution by the algorithm of table 2.1(b). Sub- 
routine MINC is called to perform this task (see fig. A. 4), and its output is shown 
in figure C7. Convergence properties are indicated by the behavior of 



INITIAL CONDITIONS: 


VARIABLE 

VALUE 


STDEV 

X2 , 

M/S2 

-0. 1205E+00 


0. OOOOE+OO 

Y2 , 

M/S2 

-0 . 8895E-01 


0 . 0000E+00 

H2 , 

M/S2 

0 . 2402E-01 


0. OOOOE+OO 

WX , 

M/S 

0 . 5003E+00 


0. 0000E+00 

WY , 

M/S 

0. 1078E+01 


0. OOOOE+OO 

WH , 

M/S 

0 . 8971E+00 


0. OOOOE+OO 

XI , 

M/S 

-0 . 2106E+01 

+ » - 

0. OOOOE+OO 

Y1 , 

M/S 

0 . 9074E+02 

+ . - 

0. OOOOE+OO 

HI , 

M/S 

0 . 3723E-01 


0. OOOOE+OO 

PHI, 

DEG 

0 . 5849E+00 

+ . ~ 

0. OOOOE+OO 

THT , 

DEG 

0. 1443E+01 


0. OOOOE+OO 

PSI, 

DEG 

0 . 9175E+02 


0. OOOOE+OO 

X , 

NM 

-0 . 9120E+00 


0. OOOOE+OO 

Y , 

NM 

-0 . 5467E+00 

+ . ~ 

0. OOOOE+OO 

H , 

M 

0 . 9999E+03 


0 . OOOOE+OO 

PHI, 

D/S 

-0 . 3503E+00 


0 . OOOOE+OO 

TH1 , 

D/S 

-0 . 5845E-02 


0 . OOOOE+OO 

PSI, 

D/S 

-0 . 3616E-01 


0. OOOOE+OO 

BIAS 

ERRORS 

: 



VARIABLE 

VALUE 


STDEV 

AX , 

G'S 

0 . 0000E+00 

+ , - 

0. OOOOE+OO 

AY , 

G'S 

0 . 0000E+00 


0 . OOOOE+OO 

A2 , 

G'S 

0 . OOOOE+OO 


0 . OOOOE+OO 

P , 

D/S 

0 . 0000E+00 

+ . - 

0 . OOOOE+OO 

9 . 

D/S 

0 . 0000E+00 

+ , - 

0. OOOOE+OO 

R , 

D/S 

0 . 0000E+00 


0. OOOOE+OO 

SCALE FACTORS: 




VARIABLE 

VALUE 


STDEV 

P 

. D/S 

0. 1000E+01 


0 . OOOOE+OO 

Q 

, D/S 

0. 1000E+01 

+ , - 

0 . OOOOE+OO 

R 

, D/S 

0. 1000E+01 


0. OOOOE+OO 


Figure C5.- Concluded. 


value of the performance measure), and GRAD (its gradient), printed at each itera- 
tion. The other parameters displayed, GAIN and DETM, are, respectively, the 
multiplier of the step and the determinant of the information matrix. Following the 
solution by MINC, control is returned to SMACK, which tabulates the solution as 
shown in figure C8. Listed there are the statistics (mean and rms values) of the 
residuals (response, or "fit," errors) and forcing-function time histories. A 
"good" solution is judged not only by robust convergence, but also by how closely 
rms values in the STDEV columns match corresponding values in the columns headed 


1 1 1 


(a) DISPLAY OF STATE MODEL ARRAYS AT T=0 : 


ARRAY: XD 


-0. 1205E+00 
-0.2333E-02 
0. 1021E-01 


0 . 2169E-01 
0 . 0000E+00 


0 . 2169E-01 
0 . 0000E+00 
-0 . 2333E-02 


-0 . 8895E-01 
-0 . 4187E-03 
0 . 2518E-01 


-0 . 5074E-02 
0 . 0000E+00 


-0 . 5074E-02 
0 . 0000E+00 
-0 . 4187E-03 


0 . 2402E-01 
0. 1182E-02 
0. 1601E+01 


-0. 1343E-02 
0 . 0000E+00 


0 . 0000E+00 
0 . 5003E+00 
-0 . 2106E+01 
-0. 1689E+04 

ARRAY: WD 


ARRAY: FD 

0 . 0000E+00 
-0. 1469E+00 
-0 . 1205E+00 
-0 . 2106E+01 

ARRAY: FX 

0. 1000E+01 
-0 . 1000E+01 
-0. 1000E+01 
0. 1210E-04 
0 . 2977E-04 


0 . 0000E+00 
0. 1078E+01 
0 . 9074E+02 
-0. 1012E+04 


0 . 0000E+00 
-0 . 5226E-01 
-0 . 8895E-01 
0 . 9074E+02 


0. 1000E+01 
-0 . 2571E-03 
0 . 2362E-02 
-0. 1186E-02 


0 . 0000E+00 
0 . 8971E+00 
0 . 3723E-01 
0 . 9999E+03 


0 . 0000E+00 
0 . 6128E-01 
0 . 2402E-01 


0 . 0000E+00 
0 . 6128E-01 
0 . 2402E-01 
0 . 3723E-01 


0. 1000E+01 
-0. 9999E+00 
0. 1045E-06 
-0. 1055E-04 


0 . 1000E+01 
-0. 1021E-01 
0 . 4066E-03 
-0 . 1182E-02 


0. 1000E+01 
-0. 2519E-01 
0 . 4151E-05 
-0.4188E-03 


-0. 1343E-02 
0 . 0000E+00 
0. 1182E-02 


0. 1000E+01 
0. 1021E-01 
-0 . 2987E-04 
0. 1182E-02 


0 . 1000E+01 


ARRAY: FW 

0.1000E+01 0.1000E+01 0.1000E+01 0.1000E+01 0.1000E+01 


0.0000E+00 0 . 0000E+00 
-0 . 1469E+00 -0 . 5226E-01 
-0 . 1205E+00 -0 . 8895E-01 


Figure C6.- Display of arrays for flight-test problem. 


(a) State model arrays at 


(R**0.5) and (Q**0.5). For this problem, the rms values match closely. Notice that 
the mean of the response error for the accelerometer fits is quite small, which 
should be the case when instrument bias errors are identified. Notice also that 
state-variables (PHI, TH1 , PS1 ) are included in the forcing-function list. These 
were computed from the measured (P, Q, R) (adjusted for bias-error and scale-factor 
parameters). Finally, the parameters identified in the solution are listed. Their 
standard deviations are found from the diagonal elements of the information matrix 
inverse, which are estimates of the Cramer-Rao lower bound. The last step m this 
solution is to return control to subroutine MODL, where a simulation summary is 
computed and printed, as shown in figure C9. A comparison of response error and 
measurement noise can be made for each variable fit to a measurement record. The 
error between the "true" solution and the estimated solution is given for each 
variable not measured. 



(b) DISPLAY OF OUTPUT MODEL ARRAYS AT T=DT: 
ARRAY: HD 


0 . 7874E-02 
0 . 2170E+04 
0 . 8967E+02 
-0 . 4904E-02 
0 . 0000E+00 

ARRAY: SD 

0. 8962E+02 
0. 1525E+00 
0. 0000E+00 
0. 0000E+00 
0. 0000E+00 

ARRAY: HX 

0. 1000E+01 
0.4608E+00 
0.2877E-01 
0. 1115E-01 
0. 9122E-04 
0. 0000E+00 
0 . 3480E-01 
-0. 9995E+00 
0 . OOOOE+OO 
-0 . 2428E-01 


0 . 2476E-01 
-0 . 2643E+01 
0. 3479E-01 
0. 1831E-04 
0. 0000E+00 


-0 . 2066E+00 
0 . 2428E-01 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


0. 1000E+01 
0 . 2485E-03 
-0 . 9995E+00 
0 . 3219E-03 
-0. 1090E-03 
0. 1421E-13 
0 . 7794E-02 
-0 . 3116E-01 
0 . 9829E+01 
0. 1525E+00 


0. 1602E+01 
0 . 4789E+00 
-0 . 2306E-02 
-0 . 8545E-03 
0. OOOOE+OO 


0. 3119E+01 
-0 . 9829E+01 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


0. 1000E+01 
-0 . 4559E-03 
0. 1001E-01 
0 . 8720E-04 
-0. 1115E-01 
0. OOOOE+OO 
-0. 1001E+01 
-0 . 7872E-02 
0. 1016E+00 
0. 1802E-02 


-0. 1691E+04 
0. 1085E+01 
0. 1525E+00 
0. OOOOE+OO 
0. OOOOE+OO 


0 . 8962E+02 
0 . 8962E+02 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


0. 1000E+01 
0. 1864E-03 
-0 . 9122E-04 
-0 . 2877E-01 
-0. 1115E-01 
0 . 2303E-02 
-0 . 3134E-01 
0 • 7094E-02 
-0 . 9829E+01 
0. 1000E+01 


-0. 9218E+03 
-0. 1903E+01 
0 . 2428E-01 
0. OOOOE+OO 
0. OOOOE+OO 


-0 . 2066E+00 
-0. 2066E+00 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


-0.7793E+00 
0. 1016E-03 
0. 1090E-03 
0 . 9995E+00 
-0 . 3219E-03 
0. 1000E+01 
0 . 9992E+00 
0 . 2499E-01 
0. 1201E-02 
0. 1000E+01 


0 . 9999E+03 
0 . 9584E+00 
-0. 9829E+01 
0. OOOOE+OO 
0. OOOOE+OO 


0. 3119E+01 
0. 3119E+01 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


-0 . 4248E+00 
0 . 4090E-03 
0. 1115E-01 
-0. 1001E-01 
-0 . 8720E-04 
0 . 7895E-02 
0 . 2476E-01 
-0 . 9997E+00 
0. 9090E-01 
0. 1000E+01 


Figure C6.- Concluded, (b) Measurement model arrays at t = DT. 


fin^l comment should be made here: the fact that rate— gyro measurements 

have been chosen to generate forcing functions in this problem, instead of estimat- 
ing (PH2, TH2, PS2) or (DL, DM, DN) and fitting estimates of (P, Q, R) to the rate- 
gyro records, probably is not important because the noise level on those records is 
quite small. The results would differ, however, if the measurements of (P, Q, R) 
were relatively noisy. The user can readily explore this topic by replacing the VAR 
statements for P, Q, and R with 

P 111 

Q 111 

R 111 

PH2 2 1 

TH2 2 1 

PS2 2 1 

and obtain another solution, and then repeat the pair of experiments with a higher 
rate-gyro noise level. A similar argument would hold for using accelerometer mea- 
surements (AX, AY, AZ) to generate forcing functions (X2, Y2, H2). In general, how- 
ever, it is preferable to fit all measurements, and estimate all forcing functions. 


1 0.001 

1 0.001 

1 0.001 
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PROBLEM SOLUTION: 

COST- 0 . 7578E+05 
GRAD- 0 . 5463E+08 

* * * XT- 1 * * * 

COST- 0 . 6617E+03 
GRAD- 0 . 2696E+06 

* * * XT- 2 * * * 

COST- 0 . 1637E+02 
GRAD- 0 . 1423E+05 

*** IT- 3 *** 

COST- 0 . 1612E+02 
GRAD- 0 . 5953E+02 

* * * it- 4 * * * 

COST- 0 . 1612E+02 
GRAD- 0 . 1163E+00 

* * * IT— 5 * * * 

COST- 0 . 1612E+02 
GRAD- 0 . 5823E-04 


GAIN- 0. OOOOE+OO 
DETM- 0 . 4058E-11 


GAIN- 0 . 1000E+01 
DETM- 0.4058E-11 


GAIN- 0 . 1000E+01 
DETM- 0.4056E-11 


GAIN- 0 . 1000E+01 
DETM- 0.4057E-11 


GAIN- 0 . 1000E+01 
DETM- 0 . 4057E-11 


GAIN- 0 . 1000E+01 
DETM- 0 . 4057E-11 


gure C7.- Iterative solution of flight-test problem. 



SOLUTION TABULATION: CPUTIME- 1.2201 SEC 
RESPONSE ERRORS AND WEIGHTS: 


VARIABLE 

MEAN 

STDEV 

(R* *0 . 5) 

PHI, 

DEG 

-0. 1118E-03 

+,- 0 . 4982E-01 

0 . 5000E-01 

THT , 

DEG 

0. 3307E-03 

+,- 0 . 4845E-01 

0 . 5000E-01 

PSI, 

DEG 

-0 . 8621E-03 

+,- 0 . 4942E-01 

0 . 5000E-01 

H , 

M 

-0. 1540E-03 

+,- 0 . 8965E-01 

0. 1000E+00 

RNG, 

NM 

-0 . 2051E-03 

+,- 0 . 4948E-02 

0 . 5000E-02 

BRG, 

DEG 

0 . 8718E-03 

+,- 0 . 4913E-01 

0 . 5000E-01 

ELV, 

DEG 

0. 1558E-02 

+ ,- 0 . 5097E-01 

0 . 5000E-01 

VT , 

XT 

0. 1254E-02 

0 . 3687E-01 

0 . 1000E+00 

AV , 

DEG 

-0. 5596E-03 

+,- 0 . 3336E-01 

0 . 5000E-01 

BV , 

DEG 

-0. 1288E-02 

+,- 0 . 2725E-01 

0 . 5000E-01 

AX , 

G'S 

-0 . 7338E-14 

+,- 0 . 7128E-02 

0. 1000E-01 

AY , 

G'S 

0 . 6943E-15 

+,- 0 . 7489E-02 

0. 1000E-01 

AZ , 

G'S 

-0. 1932E-14 

+,- 0 . 9844E-02 

0. 1000E-01 

FORCING FUNCTIONS AND WEIGHTS: 


VARIABLE 

MEAN 

STDEV 

(Q**0.5) 

DX , 

M/S3 

0. 5323E-03 

+,- 0 . 1735E+00 

0. 1643E+00 

DY , 

M/S3 

-0. 1282E-03 

+,- 0 . 9976E-01 

0 . 9321E-01 

DH , 

M/S3 

-0 . 3163E-03 

+,- 0 . 9856E-02 

0 . 9736E-02 

GX , 

M/S2 

0. 1089E-02 

+,- 0.7811E-01 

0 . 8128E-01 

GY , 

M/S2 

-0. 3480E-02 

+,- 0 . 1017E+00 

0 . 9604E-01 

GH , 

M/S2 

0. 1208E-03 

+ ,- 0 . 5334E-01 

0 . 4936E-01 

PHI, 

D/S 

-0. 1099E-03 

+,- 0 . 7382E+00 

0.0000E+00 

TH1 , 

D/S 

0 . 2510E-01 

+,- 0 . 1311E+00 

0.0000E+00 

PSI, 

D/S 

-0 . 2031E+01 

+,- 0 . 9960E+00 

0.0000E+00 


FILTER 

0. 1500E+00 
0. 1000E+00 
0. 1000E+00 
0. 1000E+00 
0. 1000E+00 
0. 1000E+00 
0. 1000E+00 
0. 1000E+00 
0. 1000E+00 
0. 1000E+00 
0. 1000E+00 
0. 1000E+00 
0. 1000E+00 


Figure C8.- Solution tabulation for flight-test problem. 


INITIAL CONDITIONS: 


VARIABLE VALUE STDEV 


X2 , 

M/S2 

0 . 8228E-02 


0 . 1789E+00 

Y2 . 

M/S2 

-0 . 6534E-01 


0. 1118E+00 

H2 , 

M/S2 

0. 1010E-01 

+ , “ 

0. 1851E-01 

WX , 

M/S 

0. 1684E+01 

+ >“ 

0. 2110E+00 

WY , 

M/S 

-0 . 4291E+01 


0 . 2040E+00 

WH , 

M/S 

0 . 9498E+00 


0 . 7920E-01 

XI , 

M/S 

-0 . 3882E+00 

+ » ~~ 

0 . 2849E+00 

Y1 , 

M/S 

0 . 8535E+02 


0 . 2345E+00 

HI , 

M/S 

-0 . 5937E-02 


0 . 5164E-01 

PHI, 

DEG 

0 . 2150E-01 


0. 1561E-01 

THT , 

DEG 

0. 1468E+01 


0. 1443E-01 

PSI, 

DEG 

0 . 9130E+02 

+ * “ 

0 . 1908E-01 

X , 

NM 

-0 . 9101E+00 


0 . 8360E-03 

Y , 

NM 

-0 . 5407E+00 


0.6059E-03 

H , 

M 

0. 1000E+04 


0. 1168E+00 

PHI, 

D/S 

-0. 1335E+00 


0 . OOOOE+OO 

TH1 , 

D/S 

-0 . 2478E-01 

+ » “ 

0 . OOOOE+OO 

PSI, 

D/S 

0 . 6582E-01 


0. OOOOE+OO 

BIAS 

ERRORS 

: 



VARIABLE 

VALUE 


STDEV 

AX , 

G'S 

0 . 1281E-03 


0 . 1111E-02 

AY , 

G'S 

0. 1310E-02 


0. 1118E-02 

AZ , 

G'S 

0 . 3850E-03 

+ , - 

0 . 1059E-02 

P , 

D/S 

-0. 1228E-03 


0 . 3238E-03 

Q . 

D/S 

0 . 1412E-02 


0. 1332E-02 

R , 

D/S 

0 . 2100E-02 

+ , “ 

0 . 2120E-02 


SCALE FACTORS: 


VARIABLE VALUE STDEV 

P , D/S 0 . 1000E+01 +,- 0 . 7339E-03 
Q , D/S 0 . 9984E+00 +,- 0.1526E-02 
R , D/S 0 . 1001E+01 +,- 0 . 9246E-03 


Figure C8.- Concluded. 
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SIMULATION RESULTS: 


VARIABLE 

PHI, DEG 
THT, DEG 
PSI, DEG 
X , NM 
Y , NM 
H , M 
RNG, NM 
BRG, DEG 
ELV, DEG 
WXY, KT 
WHD, DEG 
VWD, M/S 
VT , KT 
AV , DEG 
BV , DEG 
AX , G'S 
AY , G'S 
AZ , G'S 


ESTIMATION ERROR 

MEAN STDEV 

-0.11 18E-03 +,- 0 . 4982E-01 
0 . 3307E-03 +,- 0 . 4845E-01 
-0.8621E-03 + ,- 0.4942E-01 
0 . 1260E-03 +,- 0 . 4016E-03 
-0 . 2051E-04 +,- 0 . 3354E-03 
-0 . 1540E-03 +,- 0 . 8965E-01 
-0 . 2051E-03 +,- 0 . 4948E-02 
0 . 8718E-03 +,- 0 . 4913E-01 
0 . 1558E-02 +,- 0 . 5097E-01 
0 . 7110E-01 + ,- 0 . 2348E+00 
0 . 5131E-01 + ,- 0 . 2687E+01 
0 . 9985E-02 +,- 0.3355E-01 
0.1254E-02 +,- 0 . 3687E-01 
-0 . 5596E-03 +,- 0.3336E-01 
-0 . 1288E-02 +,- 0 . 2725E-01 
-0 . 7338E-14 +,- 0 . 7128E-02 
0 . 6943E-15 +,- 0 . 7489E-02 
-0 . 1932E-14 +,- 0 . 9844E-02 


Figure C9-- Summary of flight-test 


MEASUREMENT NOISE 

MEAN STDEV 

-0 . 3289E-02 +,- 0.5000E-01 
-0 . 2518E-02 +,- 0 . 5000E-01 
-0 . 1502E-02 +,- 0 . 5000E-01 
0.0000E+00 +,- O.OOOOE+OO 
0 . OOOOE+OO +,- 0.0000E+00 
0 . 1099E-01 +,- 0 . 1000E+00 
-0 . 1634E-03 +,- 0 . 5000E-02 
-0 . 1458E-02 +,- 0 . 5000E-01 
0 . 6007E-03 +,- 0 . 5000E-01 
0.0000E+00 +,- 0.0000E+00 
O.OOOOE+OO +,- 0.0000E+00 
0.0000E+00 +,- O.OOOOE+OO 
0 . 4505E-02 +,- 0 . 1000E+00 
-0 . 1008E-01 +,- 0 . 5000E-01 
0 . 5341E-04 +,- 0 . 5000E-01 
-0.1700E-03 +,- 0 . 1000E-01 
0 . 6972E-03 +,- 0.1000E-01 
0 . 4576E-03 +,- 0.1000E-01 


problem results. 
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APPENDIX D 


ACCIDENT PROBLEMS 


This appendix describes how SMACK can be applied for accident analysis, and 
illustrates the special output format useful for display of the results. In 
appendix A (see fig. A2), it was indicated that this option will be chosen when the 
flag ISOLU= 1 (obtained by setting K=1 in the solution description). The purpose of 
this option, which is produced by subroutine SOLU, is to provide printed time his- 
tories of trajectory information, including groundtrack, lift and drag, aircraft 
attitude, indicated airspeed, and angle of attack. Two accident problems will be 
given here, each using the simulated trajectory described in chapter 4 and 
appendix C to provide measurement records for SMACK analysis. 


Problem 1 

For this first problem, the trajectory data are obtained from two separate 
tracking stations. There are no records of onboard measurements available, except 
for barometric altitude (from the aircraft transponder), a record that is usua y 
included with ATC tracking data. Other information available to the analyst con- 
sists of air temperature and winds (from balloon measurements made as close as 
possible to the time and location of the accident), magnetic variation, and aircraft 
performance parameters. The latter include zero-lift angle of attack, wing loading, 
and the derivative of the normal-force coefficient, which are necessary for estima- 
tion of angle of attack. These parameters are normally defined when the user pre- 
pares the DATA subroutine, as described in appendix E. For the SMACK internal 
simulation, the performance parameters correspond to those of a large commercia 

airliner . 

A coding list for this problem is shown in figure D1(a). Notice that the first 
directive specifies the accident solution format (K=1 in the MKS statement). The 
locations of the tracking stations with respect to a user-selected origin are input 
next followed by the record information. Notice also that no attitude vana es 
have* been included in the coding list. Since there are no onboard measurements 
available (except altitude), the solution will be more efficient if only the point- 
mass part of the rigid-body model is excited. The summary of the coding list shown 
in figure D1(b) illustrates the form of the state and measurement models for this 
problem. In this special case, the Euler angles are determined following the fits 
to the radar data, by an algebraic method (ref. 47). This method, coded into sub- 
routine RADR, assumes zero side force and sideslip angle, and uses the performance- 
parameter derived estimate of angle of attack and the smoothed accident trajectory 
to calculate the aircraft attitude along the flightpath. The results of the angle 
calculations are included with the accident output format shown in figure D1(c). 
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(a) ACCIDENT ANALYSIS-MULTIPLE RADAR TRACK 


MKS 5 

REC 1 

RAD 

RNG 1 

BRG 1 

BRG 4 

H 1 

H 4 

WXY 1 

WHD 1 

VWD 1 

RDA 

RNA 1 

BRA 1 

BRA 4 

ELA 1 

ELA 4 

DX 2 

DY 2 

DH 2 

END 


1 1 
90 1 1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

1 


1.0 

0. 

0.005 

0.001 

0.4 

0.2 

0.075 


1. 

0.005 

0.005 

0.2 

0.005 

0.2 


0 . 100 . 


1 . 200 . 


(b) SUMMARY OF CODING LIST: 

NF - 9 (DIFFERENTIAL EQUATIONS IN STATE MODEL) 

NX - 9 (ESTIMATED VARIABLES IN STATE MODEL) 

NW - 3 (ESTIMATED FORCING FUNCTIONS) 

NH - 8 (ESTIMATED VARIABLES IN OUTPUT MODEL) 

NZ - 9 (AVAILABLE MEASUREMENTS) 

NV - 6 (OUTPUT ESTIMATES FIT TO MEASUREMENTS) 

NZH= 3 (OUTPUTS MEASURED, NOT ESTIMATED) 

NHZ- 11 (NH + NZH) 

NZX- 3 (STATES MEASURED, NOT ESTIMATED) 

NXZ- 12 (NX + NZX) 

NXO- 9 (STATE- VARIABLE INITIAL CONDITIONS) 

NWB- 0 (FORCING-FUNCTION MEAN VALUES) 

NVB- 0 (INSTRUMENT BIAS ERRORS) 

NSF- 0 (INSTRUMENT SCALE FACTORS) 

NT - 9 (NXO + NWB + NVB + NSF) 


Figure D1.- Accident problem 1 (multiple radar track), (a) Coding list, (b) list 

summary. 


Problem 2 

A second problem illustrates an accident solution for a measurement set that 
includes not only tracking data (from a single station) and winds, but also onboard- 
instrument data records typically recovered from a metal-foil crash recorder. Such 
records contain time histories of airspeed, pressure altitude, magnetic heading and 
vertical acceleration. For this problem, it is assumed that the heading measurement 
comes from a two-degree-of-freedom directional gyro of the type described in 
appendix B. A coding list for the measurement set is shown in figure D2(a). 
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(c) ACCIDENT ANALYSIS RESULTS FOR PROBLEM 1 


OINT 




GROUND 

TRACK 

NO 


ALTITUDE 

SPEED 

ANGLE 


MIN 

SEC 

FT 

knots 

DEG 

1 

0 

1 00 

3280 

164 9 

90 13 

2 

0 

2 00 

3280 

165.1 

90 07 

3 

0 

3 00 

3280 

165.1 

09.95 

4 

0 

4 00 

3281 

165.0 

89.78 

5 

0 

5 00 

3281 

164.9 

89 . 58 

6 

0 

6 00 

3281 

1 64.6 

89 30 

7 

0 

7 00 

3281 

1 64, 1 

88 89 

8 

0 

8 00 

3282 

163 6 

88 32 

9 

0 

9 00 

32B2 

163 2 

87 . 60 

10 

0 

10 00 

3283 

162 6 

86 72 

1 1 

0 

1 1 00 

3284 

162.0 

85 64 

12 

0 

12 00 

3285 

161.1 

84 36 

13 

0 

13 00 

3207. 

160.1 

82 91 

1 4 

0 

14 00 

3288 

159.3 

Bi 32 

15 

0 

15 00 

3291 . 

156 3 

79 62 

16 

0 

16 00 

3294 

157,3 

77.77 

17 

0 

17 00 

3297. 

156 5 

75 69 

18 

0 

18 00 

3301 

155 7 

73 43 

19 

0 

19 00 

3305 

155 1 

70 98 

20 

0 

20 00 

3310 

154 5 

68.42 

21 

0 

21 00 

3316 

154 0 

65.77 

22 

0 

22 00 

3322. 

153.7 

63 02 

23 

0 

23 00 

3329. 

153 6 

60 27 

24 

0 

24 00 

3337. 

153.5 

57 49 

25 

0 

25 00 

3345 

153 5 

54 68 

26 

0 

26 00 

3354 

1538 

51 83 

27 

0 

27 00 

3364 

154.1 

48 98 

28 

0 

28 00 

3374 

154 3 

46.15 

29 

0 

29.00 

3385 

154 4 

43.38 

30 

0 

30 00 

3396 

154 6 

40 62 

31 

0 

31 00 

3407. 

154 7 

37 91 

32 

0 

32 00 

3420. 

154 5 

35.25 

33 

0 

33 00 

3432. 

154 3 

32 61 

34 

0 

34 00 

3445 

153.9 

30 00 

35 

0 

35 00 

3459 

153 4 

27 38 

36 

0 

36 00 

3472 

152 8 

24.75 

37 

0 

37 00 

3486 

152 0 

22.11 

38 

0 

38 00 

3500 

151.0 

19 45 

39 

0 

39.00 

3515 

150.0 

16 78 

40 

0 

40 00 

3529 

149 0 

14.10 

41 

0 

41 00 

3543 

148.0 

11.39 

42 

0 

42 00 

3558 

147.1 

8.61 

43 

0 

43 00 

3572 

146 5 

5.77 

44 

0 

44 00 

3587. 

146 1 

2 90 

45 

0 

45 00 

3602. 

145 8 

0 02 

46 

0 

46 00 

3616 

145 9 

357.12 

47 

0 

47.00 

3631 . 

146 5 

354,24 

48 

0 

48 00 

3645. 

147 . 2 

351 40 

49 

0 

49 00 

3660 

148 1 

348.61 

50 

0 

50 00 

3674 

148.9 

345.86 

51 

0 

51 00 

3689 

150.0 

343.16 

52 

0 

52 00 

3703. 

151.1 

340 50 

53 

0 

53 00 

3718 

151.9 

337 87 

54 

0 

54 00 

3732- 

152.7 

335 26 

55 

0 

55 00 

3745 

153 4 

332 66 

56 

0 

56 00 

3759. 

154,0 

330 06 

57 

0 

57 00 

3772 . 

154 3 

327.41 

58 

0 

58 00 

3785 

154 5 

324 74 

59 

0 

59 00 

3798 

154 7 

322 06 

60 


0 . 00 

3810 

154 5 

319 33 

61 


1 .00 

3822 

154 4 

316.61 

62 


2 00 

3833 

154.2 

313.85 

63 


3 00 

3844 

154.1 

31 1 06 

64 


4 00 

3854 . 

153.7 

308.17 

65 


5 00 

3863. 

153.6 

305.33 

66 


6 00 

3872. 

153.5 

302 . 52 

67 


7 00 

3881 . 

153 6 

299 75 

68 


8 00 

3888 

153.7 

296.96 

69 


9 00 

3895 

154 0 

294 24 

70 


1 0 , 00 

3902. 

154.5 

291 60 

71 


1 1 00 

3907. 

155 0 

289 02 

72 


12.00 

3912. 

155 7 

286.59 

73 


13 00 

3917 

156.6 

284 38 

74 


14 00 

3921 . 

157 4 

282.31 

75 


15.00 

3924. 

158.3 

280 42 

76 


16 00 

3927. 

159 3 

278 68 

77 


17 00 

3929. 

160 2 

277 07 

78 


18 00 

3931 . 

161 1 

275 62 

79 


19.00 

3933 

161 8 

274.36 

80 


20 00 

3934. 

162 6 

273.30 

81 


21.00 

3935. 

163.3 

272.42 

82 


22.00 

3936 

163.8 

271 .74 

83 


23 00 

3936. 

164.3 

271 ,21 

84 


24.00 

3937. 

164.7 

270.79 

85 

1 

25 00 

3937. 

165 0 

270 46 

86 

1 

26 00 

3937. 

165,1 

270.25 

87 


27.00 

3937 

165.1 

270 09 

80 


28.00 

3937 

165.2 

269 98 

89 


29 00 

3937 

165,1 

269.86 

90 


30 00 

3937 

165 1 

269 75 


VERT 

FLIGHT 

FORCES 


VEL 

PATH 

LIFT 

T-D 

ROIL 

FPS 

DEG 

G ' S 

G'S 

DEG 

0 02 

0,00 

1 00 

-0 01 

-0 56 

0 07 

0.01 

1 00 

-0 01 

-0 99 

0.12 

0 02 

1 00 

-0.02 

-1.49 

0 18 

0 04 

1 00 

-0.02 

-1 .67 

0 24 

0 05 

1 00 

-0.03 

-2 46 

0 32 

0 07 

1 00 

-0.03 

-3.50 

0 42 

0 09 

1 .01 

-0.04 

-4 78 

0.54 

0.11 

1 .01 

-0 04 

-6 06 

0. 70 

0.15 

1 01 

-0 04 

-7 45 

0.92 

0.19 

1 02 

-0 05 

-8 95 

1 . 18 

0 25 

1 03 

-0.06 

-10.53 

1 .51 

0 32 

1 03 

-0 07 

-11.86 

1 .89 

0 40 

1 04 

-0.06 

-12.85 

232 

0 50 

1 04 

-0.07 

-13.58 

2.81 

0 60 

1 .05 

-0.07 

-14.72 

3.33 

0.72 

1 .06 

-0.06 

-16.29 

3.89 

0.84 

1 07 

-0 06 

-17.57 

4 48 

0.98 

1 08 

-0 05 

-18 82 

5. 10 

1.12 

1 08 

-0.05 

-19 . 53 

5.73 

1 .26 

1 09 

-0.04 

-20. 19 

6 37 

1 . 40 

1 09 

-0.03 

-20.78 

7 01 

1 . 55 

1 .09 

-0.02 

-20 87 

7.65 

1 69 

1 .09 

-0.02 

-21 .01 

8 29 

1 83 

1 .09 

-0.01 

-2122 

8 91 

1 . 97 

1 09 

0 00 

-21 .61 

9.51 

2. 10 

1 09 

0 01 

-21 64 

10.09 

222 

1 . 09 

0 00 

-21 .47 

10 65 

2.34 

1 09 

0 00 

-21 20 

11.19 

2.46 

1 .09 

0 01 

-2111 

11.70 

2.57 

1 08 

0 00 

-20.75 

12.18 

2 67 

1 08 

-0 01 

-20.44 

12.62 

2.77 

1 08 

-0.01 

-20.34 

13 02 

2,86 

1 .08 

-0.02 

-20.04 

13 36 

2 95 

1 07 

-0 02 

-20 08 

13 65 

3.02 

1 .07 

-0 02 

-20 14 

13.89 

3.08 

1 .07 

-0.04 

-20 1 1 

14 09 

3.14 

1 07 

-0.04 

-20 22 

14.24 

3 20 

1 07 

-0 05 

-20. 13 

14.36 

3.25 

1 . 07 

-0.04 

-20.09 

1 4 44 

3.29 

1 07 

-0.05 

-20.22 

14 50 

3.32 

1 .07 

-0.04 

-20.57 

14 54 

3 35 

1 07 

-0.02 

-20.97 

14 57 

3.37 

1 07 

-0 02 

-21.12 

14 60 

3 39 

1 .07 

-0.01 

-21.20 

14 61 

3.40 

1 07 

0 01 

-21.33 

14 61 

3 . 40 

1 .07 

0 03 

-21 . 32 

14.60 

3.38 

1 07 

0 05 

-21.14 

1 4 57 

3 36 

1 07 

0 05 

-20.98 

14.52 

3 33 

1 07 

0 05 

-20.77 

1 4 45 

3.29 

1 06 

0 06 

-20.60 

14 35 

3.25 

1 06 

0 06 

-20.45 

14.23 

3 19 

1 06 

0 05 

-20.39 

14.06 

3 14 

1 .06 

0 05 

-20.31 

13.86 

3.08 

1 .06 

0 04 

-20.37 

13.61 

3.01 

1 06 

0 03 

-20.43 

13.32 

2 93 

1 06 

0 02 

-20.81 

12.99 

2 86 

1 06 

0.01 

-21 06 

12 61 

2.77 

1 06 

0 01 

-21.15 

12 18 

2 67 

1 .06 

-0 01 

-21.46 

11.72 

2.57 

1 06 

-0 01 

-21 . 44 

11.22 

2 47 

1 . 06 

-0 02 

-21 69 

10.70 

2.35 

1 . 06 

-0 02 

-21 .92 

10.14 

2.23 

1 06 

-0 03 

-22.55 

9 55 

2.11 

1 .06 

-0.02 

-22.27 

8 94 

1 .98 

1 06 

-0 02 

-22.03 

8 31 

1 84 

1 .05 

-0.01 

-21 . 79 

7.66 

1 . 69 

1 .06 

-0 01 

-21 .92 

7.01 

1 .55 

1 .05 

0 00 

-21.48 

6.35 

1 .40 

1 .05 

0 01 

-21 01 

5-71 

1 .25 

1 .05 

0 01 

-20 60 

5.07 

1.11 

1 .04 

0.02 

-19 52 

4 . 45 

0 97 

1 .03 

0.03 

-17.98 

3 86 

0,84 

1 03 

0.03 

-16 95 

3 31 

0.71 

1 .02 

0.04 

-15.69 

2.79 

0.60 

1 02 

0.04 

-14.44 

2.32 

0 50 

1 .02 

0.04 

-13.49 

1 91 

0 . 40 

101 

0.04 

-12.19 

1 54 

0.33 

1 01 

0 03 

-10 70 

1 .23 

0.26 

1 .00 

0 03 

-9 08 

0.97 

0.20 

1 . 00 

0.03 

-7 50 

0.76 

0. 16 

1 . 00 

0.02 

-5 83 

0.58 

0. 12 

1 00 

0 01 

-4 55 

0.44 

0.09 

1 00 

0 01 

-3.62 

0.31 

0.06 

1 00 

0 01 

-2.85 

0.21 

0.04 

1 .00 

-0.01 

-1.86 

0.11 

0,02 

1 .00 

-0.01 

-1 .35 

0.03 

0 .01 

1 .00 

-0.01 

-1 .02 

-0.05 

-0.01 

1 . 00 

-0.01 

-0.97 

-0.14 

-0.03 

1 .00 

-0.01 

-0.97 

-0.22 

-0 04 

1 00 

-0 01 

-0.97 


ancles airspeed 


PITCH 

HEADING 

TRUE 

IND 

ALPHA 

DEG 

DEG MAG 

KNOTS 

KNOTS 

DEG 

1 47 

91 43 

173 8 

165 8 

2 11 

1 46 

91 34 

173 9 

165 9 

2 10 

1 47 

91 . 19 

173.9 

165 8 

2 10 

1 49 

90 99 

173.8 

165.7 

2 12 

1 51 

90.73 

173 6 

165 6 

2 14 

1 56 

90.37 

173 3 

165 . 2 

2 19 

1 63 

89 87 

172 8 

164.7 

2.27 

1 72 

89 .21 

172.2 

164.2 

2 37 

1 82 

88 38 

171 6 

163 6 

2 47 

1 95 

87 38 

1 70 9 

162 9 

2 60 

2.11 

86 18 

170.1 

162.1 

2 76 

2.31 

84.77 

169.0 

161.1 

2 95 

2 52 

83 20 

167.8 

160 0 

3.14 

2 73 

81 . 51 

166 8 

159 0 

3.31 

2 98 

79 68 

165 6 

157 9 

3 53 

3 26 

77 64 

164 4 

156.7 

3 78 

3 . 52 

75 40 

163.3 

155.7 

4 01 

3 79 

72 96 

162 3 

1 54 7 

4 24 

4 02 

70 40 

161 4 

153 8 

4 42 

4 26 

67 73 

160 5 

153 0 

4 60 

4 48 

64 96 

159-7 

152 2 

4 76 

4 65 

62.17 

1 59 2 

151 7 

4 85 

4 80 

59 36 

158.8 

151 .3 

4 92 

4 95 

56.53 

156 4 

150 9 

4 99 

5 08 

53 65 

158 2 

150.7 

5 05 

5 16 

50 . 78 

158.2 

150 6 

5 06 

523 

47.93 

158.2 

150 6 

5 04 

5 30 

45 1 1 

158.1 

150 5 

5 02 

5 37 

42 32 

158.0 

150 4 

5 02 

5.42 

39 58 

158 0 

150 3 

4 99 

5 49 

36 88 

157 8 

150 2 

4 98 

5 58 

34 19 

157.5 

149 8 

5.02 

5 . 67 

31 54 

157 0 

149 4 

5 05 

5 79 

28 88 

156 4 

148.7 

5 14 

5 92 

26.21 

155 7 

148.1 

5 24 

6 04 

23 53 

155 0 

147.3 

5.34 

6 21 

20 82 

154 0 

146 4 

5.50 

6 39 

18 09 

152 9 

145 3 

5 67 

6 . 58 

15 36 

151.7 

144 1 

5.86 

6 76 

12.60 

150 6 

143 0 

6.05 

6 96 

9.78 

149.5 

142 0 

6.27 

7 15 

6.89 

148 5 

1410 

6 49 

7 28 

3 99 

147 8 

140 3 

6 62 

7 38 

1 .09 

147.3 

139 8 

6 73 

7 . 45 

358.18 

146.8 

139.3 

6 81 

7 44 

355 32 

146 9 

139 4 

6 79 

7.34 

352 52 

147.4 

139 8 

6 69 

7.21 

349.79 

148.1 

1 40 4 

6.54 

7.07 

347 1 1 

148 8 

1411 

6 38 

6.92 

344.47 

149 6 

1418 

6 22 

6.73 

341 90 

150 6 

142 7 

6 03 

6 56 

339 . 36 

151 .5 

143 6 

5 86 

6.42 

336.83 

152 2 

144 2 

5.72 

6 29 

334.32 

152.9 

144 8 

560 

6 18 

331 .80 

153 4 

145 3 

5 52 

6. 10 

329.24 
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145 6 

5 48 

6 05 

326.63 

153 8 

145 6 
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6,01 

324 01 

153 8 

145 5 

5 49 

5 99 
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153 6 

145 3 

5 53 

6.00 

316.63 
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144 8 

5 60 

6.01 

315 89 
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5 70 

6.04 
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152 0 

143 7 

5 81 

6.09 

310 19 
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6 07 
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304 40 
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1417 

6 15 
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301 53 
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141 . 2 
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6 09 
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148 9 

140 7 

6.32 

6 04 

295.79 

148 5 

140 3 

6.35 

5 96 

293.00 

148 4 

140 2 

6 35 

5 85 

290.27 

148 4 
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6 32 

5.71 

287 67 

148 5 

140 3 

6.22 

5.52 

2B5 28 

1488 

140 5 

6 06 
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283 04 

149 3 
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5.14 
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268.44 

156.4 

147.7 

4 44 

3.70 

268.30 

156.4 

147 7 

4.44 

3 68 

268.18 

156.4 

147 6 

4 44 


Figure D1.- Concluded, (c) Solution output format. 


ORIGINAL PAGE IS 
OF POOR QUALITY 



(a) ACCIDENT ANALYSIS-RADAR+ AIRSPEED, HEADING 


MKS 5 

REC 1 

RAD 

RNG 1 

BRG 1 

BRG 4 

H 1 

H 4 

WXY 1 

WHD 1 

VWD 1 

HDG 1 

VT 1 

AX 

AY 1 

AZ 1 

AV 1 

BV 1 

PHI 
THT 

PH2 2 

TH2 2 

PS2 2 

DX 2 

DY 2 

DH 2 

END 


1 1 
90 1 1 

1 

1 

1 


1 

1 

1 

1 

1 1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
2 


1.0 

0. 

0.005 

0.001 

0.4 

0.2 

0.075 


0.05 

0.1 

0.01 

0.01 

0.02 

0.02 


0. 


100 . 


(b) SUMMARY OF CODING LIST: 


NF - 15 (DIFFERENTIAL EQUATIONS IN STATE MODEL) 
NX - 15 (ESTIMATED VARIABLES IN STATE MODEL) 

NW - 6 (ESTIMATED FORCING FUNCTIONS) 

NH - 15 (ESTIMATED VARIABLES IN OUTPUT MODEL) 

NZ - 12 (AVAILABLE MEASUREMENTS) 

NV - 9 (OUTPUT ESTIMATES FIT TO MEASUREMENTS) 

NZH- 3 (OUTPUTS MEASURED, NOT ESTIMATED) 

NHZ- 18 (NH + NZH) 

NZX- 3 (STATES MEASURED, NOT ESTIMATED) 

NXZ- 18 (NX + NZX) 

NXO- 15 (STATE -VARIABLE INITIAL CONDITIONS) 

NVB- 0 (FORCING-FUNCTION MEAN VALUES) 

NVB- 1 (INSTRUMENT BIAS ERRORS) 

NSF- 0 (INSTRUMENT SCALE FACTORS) 

NT - 16 (NXO + NVB + NVB + NSF) 


Figure D2.- Accident problem 2 (radar, airspeed, heading), (a) Coding list, 

(b) list summary. 

Because there are now body-axis measurements available, the complete rigid-body 
aircraft model can be specified in the coding list. Note, however, that the data 
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(c) ACCIDENT ANALYSIS RESULTS FOR PROBLEM 2 

point GROUND TRACK 

NO ALTITUDE SPEED ANGLE 

MIN SEC FT KNOTS DEG 

1 0 1 .00 3280 165.2 89.94 

2 0 2 00 3280 1652 90.01 

3 0 3 00 3280 1652 89 95 

4 0 4 00 3281 . 165.1 89 80 

5 0 5 00 3281 . 164 9 89 59 

6 0 6 00 3281 , 164.7 89 24 

7 0 7 00 3281 164,3 88 86 

8 0 8 00 3282. 163 7 88 31 

9 0 9 00 3282. 163 3 87 54 

10 0 10 00 3283 . 1627 86 65 

11 0 11 00 3284 162 0 85 59 

12 0 12 00 3285 161 .1 84 35 

13 0 13.00 3287. 160 1 82 93 

14 0 14 00 3289 159 3 81 36 

15 0 15 00 3291 158 3 79.62 

16 0 16 00 3294 157 .3 77.66 

17 0 17 00 3297 1 56 5 75 60 

IB 0 18 00 3301 . 155 7 73.34 

19 0 19 00 3306 155 0 70 95 

20 0 20 00 3311 154 5 68 41 

21 0 21 00 3316 154 0 65.76 

22 0 22 00 3323. 153,7 63.06 

23 0 23 00 3330 153 6 60.30 

24 0 24 00 3337. 153.5 57 48 

25 0 25 00 3345 1 53 6 54.62 

26 0 26 00 3354 153.7 51.76 

27 0 27 00 3364 1 54 0 48 90 

28 0 28 00 3374 154 2 46.12 

29 0 29 00 3384 1 54 4 43.36 

30 0 30 00 3396 1 54 6 40.58 

31 0 31 00 3407. 154 6 37 88 

32 0 32 00 3419 154 5 35.21 

33 0 33 60 3432 1 54 3 32 61 

34 0 34 00 3445 153 .8 30 02 

35 0 35 00 3459 1 53 4 27 43 

36 0 36 00 3472 152 8 24 82 

37 0 37 00 3486 151 9 22.21 

38 0 38 00 3500 151.1 19.55 

39 0 39 00 3515 1 50 0 16 88 

40 0 40 00 3529 149 0 14.18 

41 0 41 00 3543 148 0 1 1 40 

42 0 42 00 3558. 147.2 8.58 

43 0 43 00 3572. 146 5 5.73 

44 0 44.00 3587. 146. 1 2.90 

45 0 45 00 3602 1 45 8 001 

46 0 46 00 3616 146 0 357 13 

47 0 47 00 3631 . 146 5 354.24 

48 0 48 00 3645. 147.2 351 40 

49 0 49 00 3660. 148 1 348 65 

50 0 50 00 3674. 149 0 345 90 

51 0 51 00 3689 150 0 343.21 

52 0 52.00 3703. 151 . 1 340.52 

53 0 53 00 3717. 151 9 33? 87 

54 0 54 00 3731 . 152.8 335.23 

55 0 55 00 3745 153.5 332.64 

56 0 56 00 3759 154 0 330.02 

57 0 57 00 3772. 154.3 327.39 

58 0 58 00 3785 154 6 324 75 

59 0 59 00 3798. 154.6 322 05 

60 1 0 00 3810 154.5 319.34 

61 1 1 00 3822 154 4 316 62 

62 1 2 00 3B33 154 3 313 88 

63 1 3 00 3844 154.0 311 05 

64 1 4 00 3854 153.8 308 22 

65 1 5.00 3863 153.6 305.36 

66 1 6 00 3872. 153.6 302.50 

67 1 7 00 3881 153 6 299 69 

68 1 8 00 3888 153.7 296 94 

69 1 9 00 3895 154 0 294 19 

70 1 10 00 3901 154.5 291.53 

71 111 00 3907 155.0 289.01 

72 1 12.00 3912. 155.7 286.64 

73 1 13.00 3917. 156.6 284 40 

74 1 14 00 3921. 157 4 282 26 

75 1 15.00 3924 158 3 280 35 

76 1 16 00 3927 159 2 278 56 

77 1 17 00 3929 160 2 276.97 

78 1 16 00 3931 1611 275.58 

79 1 19 00 3933. 161 .8 274 40 

60 1 20 00 3934. 1 62 6 273.38 

81 1 21 00 3935 163.3 272.52 

02 1 22 00 3936 163 9 271 79 

83 1 23 00 3936 1 64 3 271 .22 

84 1 24 00 3937 164.7 270.80 

85 1 25 00 3937. 164.9 270.44 

66 1 26 00 3937. 1 65 0 270.21 

87 1 27.00 3937. 165.1 270 08 

88 1 28.00 3937. 165.2 269 99 

89 1 2900 3937. 165,2 269 91 

90 1 30 00 3936 165.2 269.96 


VERT FLIGHT FORCES 

VEL PATH LIFT T-D ROLL 

FPS DEG G ' S G’S DEG 

0.00 0 00 1.00 -001 0 63 

0 06 001 1 00 -001 -0.11 

0.11 0 02 1.00 -0 02 -0 07 

015 0 03 1 00 -0.02 -1 67 

0.21 0.04 1.00 -0.02 -2. 52 

0 28 0.06 1.00 -0.03 -3 45 

0 39 0.08 1.01 -0 04 -4.47 

0 54 0 11 1.01 -0 04 -5.57 

0 73 0. 15 1 .01 -005 -6.76 

0.95 0 20 1 02 -0.05 -8.03 

1 23 0 26 1 03 -0 06 -9 35 

1 56 033 1 03 -007 -1071 

1.93 041 1 04 -0 06 -1 2 07 

2.36 0 50 1 04 -0 07 -13 40 

2.83 0 61 1 05 -0 07 -14 69 

3.36 0 72 1 06 -0 06 -15 91 

3 92 0 85 1 07 -0 06 -17 04 

4 50 0 98 1.07 -0.05 -18 08 

5.10 1.12 1.08 -0.05 -1901 

5.71 1 .26 1 09 -0 04 -19 80 

6.33 1.40 1.09 -0 03 -20 42 

6 96 1.54 1.09 -0 02 -20.86 

7.59 1.68 1 09 -001 -21.14 

8 23 1.82 110 -001 -21 26 

8 85 1.96 110 0 00 -21.29 

9 46 2.09 1 09 0 01 -21 25 

10 05 2.21 1.09 0 00 -21.18 

1062 2.34 1 09 001 -21.09 

1117 2.45 1.09 0 01 -20 99 

11.69 2.57 1.08 0 00 -20.88 

12.18 2.67 1 .08 000 -20.75 

12 .64 2.77 1 08 -0.01 -20.61 

13.04 2.87 1.07 -0.02 -20 48 

13 39 2.95 1.07 - 0.02 -20 38 

13.69 3.03 1.07 -0 02 -20.30 

13.93 309 1 07 -0 04 -20.28 

1413 3.16 1.07 -0 04 -20 29 

1 4 29 3.21 1 07 -0 05 -20 35 

14 39 325 1.07 -0.05 -20.48 

14 46 3 29 1.07 -0 04 -20.64 

14 50 3.32 1 07 -0.04 -20.81 

14 53 3.35 1.07 -003 -2097 

14 54 3.37 1.07 -0.02 -21.12 

14 54 3.38 1 07 -001 -21 .22 

14.53 3 38 1.07 0 01 -21 .23 

14 54 3 38 1 07 0 03 -21 14 

1 4 54 3 37 107 0 04 -20-98 

14.53 3.35 1.07 0.05 -20.79 

14 51 3 32 1 07 0.05 -20 59 

14 47 3 29 1.06 0.06 -20.42 

14 39 3.25 1.06 0 06 -20.28 

14.28 3.21 1.06 0.05 -20 18 

14.12 3.15 1 06 0 05 -20 14 

13,92 3 09 1 06 0 04 -20 1 5 

13.67 3 02 1 06 0 03 -20.21 

13.38 2 95 106 002 -20.32 

1304 2.87 1.06 0.01 -20.49 

12 66 2 78 1 06 0.00 -20.74 

12.22 2 68 1 06 -0.01 -21 04 

11.73 2.58 106 -0.01 -21.41 

11,21 2.46 1 06 -0,02 -21 76 

10.65 2-34 1 06 -0 02 -22 05 

10.07 2.22 1 06 -0 02 -22 26 

9 47 2 09 1 06 -0,02 -22.39 

6.85 1 96 1 06 -0.01 -22 .41 

8.23 1 82 106 -0.02 -22.33 

7 60 1 68 1 .05 -0.01 -22 1 1 

6 97 1.54 1.05 0 00 -2 1 73 

6.35 140 1.05 001 -2119 

5.72 1.26 1 04 0 01 -20 52 

5*11 1.12 1.04 0 02 —19.71 

4.52 0 99 1 03 0 03 -18 77 

3.95 0 86 1 03 0 03 -17.70 

3 40 0 73 1-02 0 04 -16.55 

2 88 0.62 102 0 03 -15 32 

2.40 0.51 1 01 0.04 -14 02 

1.96 041 1.01 0 04 -12.66 

1 .56 0 33 1 01 0 03 -11 27 

1.23 0 26 1 00 0 03 -9.86 

0 95 0 20 1 00 003 -8.46 

0 72 0 15 1.00 0 02 -7,09 

0.54 011 1 .00 001 -5.79 

0.38 0.08 1 00 001 -4.56 

0.24 0 .05 1 00 000 -3.44 

011 0 02 1.00 -0 01 -2.44 

-0 01 0.00 1.00 -0.01 -1.57 

-0 12 -0 03 1.00 -0.01 -0.80 

-0 23 -0.05 1.00 -0.01 -0.11 

-0 34 -0.07 1.00 -0.02 053 

-0 44 -0 09 1 00 -002 1.16 


ANGLES AIRSPEED 

PITCH HEADING TRUE I NO ALPHA 
DEG DEG MAG KNOTS KNOTS DEG 

1 45 91 31 174 1 166 0 2 09 
1 44 91 .31 174 1 166 0 2-08 
1 45 91 19 174 0 165 9 2 08 
1 47 91 .02 1739 165 8 211 
1 52 9076 173.7 165-6 2 16 

1,57 90 33 173 4 1653 222 
1.63 89 82 172 9 1 64 8 2 28 
1 70 89 19 172.2 164.2 2 36 
1 81 88 35 171 .7 1 63 7 2 45 
1.95 87.34 170 9 1630 2 59 
211 86 19 170.1 162.1 2 74 
2.29 84 85 169 0 161 1 2 90 
2.50 83. 2B 167.9 1 60 1 3 09 
2.74 B1 58 166.8 159.1 3 30 
3 00 79.66 165 6 1 57 9 3.55 
3 25 77.58 164 4 1 56 7 3 76 
3 51 7536 1633 155 7 3 97 

3 76 72-96 1 62 . 3 154 7 4 16 

4 01 70 43 161 .3 153-8 4 38 
4 24 67 79 160 5 153 0 4 56 
4 46 65 03 159 7 1 52 2 4 71 
4 65 62.22 159-2 1517 4 85 
4 82 59.37 158 8 1513 4 95 

4 96 56.52 1 58 5 1510 501 

5 07 53.62 158 2 1 50 7 5 04 
5.16 50.72 158 1 1 50 6 5 06 

5.23 47.87 1 58 1 1 50 6 5 04 
5 30 45.11 158.1 150 5 5. 02 
5 37 42 32 1 58 0 1 50 4 5 02 
5 43 39 55 1 57 9 1 50 3 5 00 
5 50 36 81 157.7 1 50 1 501 

5.57 34 13 157.4 149 8 5 02 
5 66 31 50 157 0 149 4 5.06 
5 78 28 88 1 56 3 1 48 7 5 13 

5 91 26.23 155.7 148.1 523 
6.06 23.62 154 9 147.3 5 34 

6 21 20.91 1539 146.3 5 50 
6 38 18.18 152.9 1 45 3 5 66 

6.58 15.42 151 .7 144 2 5 86 
6 78 12.61 1 50 6 1 43 0 6 09 

6 97 9.75 1 49 5 142 0 6 29 
7. 13 6.87 148 5 141 0 6 47 
7.27 3.96 147.8 140.3 6 62 

7.37 1.09 147.2 139.7 6.73 
7.42 35821 146 9 139 4 6 78 

7 40 355 34 1 47 0 1 39 4 6 77 
7,33 352.55 147.4 139 8 6.67 
7 22 349 80 148.0 140 4 6.55 
7 08 347.14 148 9 1412 6 39 
6 92 344.50 149 6 141.8 6.23 
6 76 341.97 150 6 142 7 6 04 
6 60 339.40 151.5 143.6 5 88 
6.45 336. B3 152 3 144.2 5 74 
6.32 334.31 152 9 144 8 5 61 
6 20 331.77 153.4 1 45 3 5 53 

6.11 329.27 153.8 145 6 5 46 
6 05 326 67 153.9 145 6 5 44 
6 00 324 04 1 53 8 145 6 546 
5 98 321.37 153 6 145 3 5 50 

5 98 318.66 153 1 144 8 5.58 

6 00 315.91 152.6 144 4 5 68 
6 03 313. 12 152 0 143 8 5 81 
6 06 31023 151.3 143.1 5.9* 

6 10 307.29 150.6 142.4 6 09 

6.11 304 39 1 50 0 1418 6.19 

6.11 301 47 149.5 141.2 6.26 
6 07 298 57 1 48 9 1 40 7 6 32 
6 03 295.71 1 48 6 140 4 6 37 
5 95 292.90 148 4 140.2 6 36 
5.86 290.20 148.4 1 40 2 6 32 
5 72 287 62 148.5 140 3 6 24 
5 55 285.21 1 48 8 140.5 6. 12 
5.36 282 95 1 49 3 141 0 5.96 
5 16 280 80 149.9 141.5 5.79 
4 96 278 87 150 6 142.2 5.62 
4 78 277 .10 151.2 142 8 5 45 
4.60 275.47 152 0 143.5 5 30 
4 42 274 07 152 8 144 2 5 14 
4 26 272 86 153.4 144 8 4 99 

4.12 271 87 1 54 0 1 45 5 4 85 
4 00 271.02 154.7 146.1 4 73 
3 90 270 26 155.2 146 6 4 63 
3 83 269 66 155 6 146.9 4 57 
3 78 269.24 156 0 147 3 4 52 
3.75 268 92 156 2 1 47 5 4 49 
3 73 268 71 156.3 147.6 4 48 
3.71 268.57 1 56 4 147 6 4 47 
3 69 268.52 1 56 4 147.7 4 46 
3 68 268.47 156.5 147.7 4 46 
3 65 268 51 1 56 4 147.7 4 46 


Figure D2.- Concluded. 


(c) Solution output format. 


ORIGINAL PAGE IS 
OF POOR QUALITY 
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set is still not sufficient to allow the removal of the assumption of zero side 
orce and sideslip angle. For the simulated trajectory, which consists of a fully- 
coordinated turn, the assumption is justified. The analyst puts these constraints 
on the problem by declaring measurements of AY and BV as indicated in the coding 
ist, and placing zeros in the corresponding time-history arrays (see appendix E) 
The coding-list summary given in figure D2(b) shows that the proper attitude and 
position variables are excited in the state and measurement models. 

The ° Utp ^ fc f °™ afc for this P r °blem is shown in figure D2(c). One might expect 
better estimates of the roll and pitch angle time histories in this problem compared 
to the previous case because of the additional measurements. The improvement 
actually realized, however, depends on the quality of the additional measurements 
as well as on the weights assigned to the "pseudo" measurements (AY and BV) used to 
match the assumptions made in solving the first problem. A comparison of the two 

Sure°D^ ^ y .° bservin g the Emulation results for the two cases shown in 

oi? T u 3 ’ - Th largest improvement is seen to be in the estimate of roll angle 

' „ TblS ^/Hustrated by the roll-angle plots of the estimates compared with the 
true time histories shown in figure D4. 
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(a) 

VARIABLE 

X , NM 
Y , NM 
H , M 
RNG, NM 
BRG, DEG 
RNA, NM 
BRA, DEG 
ELA, DEG 
PHI, DEG 
THT , DEG 
PSI, DEG 
VT , XT 
AV , DEG 
AX , G'S 
AZ , G'S 


(b) 

VARIABLE 

PHI, DEG 
THT, DEG 
PSI, DEG 
X , NM 
Y , NM 
H , M 
RNG, NM 
BRG, DEG 
VT , XT 
AV , DEG 
BV , DEG 
AX , G'S 
AY , G'S 
AZ , G'S 
HDG, DEG 


Figure 


ESTIMATION ERROR 

MEAN STDEV 

0 . 3736E-04 +,- 0.1092E-03 
0 . 2122E-04 +,- 0 . 4867E-04 
0 . 2749E-02 +,- 0.1952E+00 
0 . 2279E-03 +,- 0.5026E-02 
■0.1131E-04 +,- 0 . 4740E-03 
0 . 5844E-03 +,- 0.4995E-02 
0 . 5934E-04 +,- 0.4028E-02 
-0 . 1537E-03 +,- 0 . 4759E-02 
0 . 3442E-01 +,- 0 . 4432E+00 
-0 . 7822E-03 +,- 0.1983E-01 
-0 . 6208E-02 + ,- 0 . 6060E-01 
0 . 9709E-02 +,- 0 . 6211E-01 
0 . 1251E-02 + ,- 0 . 2247E-01 
-0 . 8864E-04 +,- 0.3371E-02 
-0 . 5054E-04 +,- 0 . 1294E-02 


MEASUREMENT NOISE 

MEAN STDEV 

0 . 0000E+00 +,- 0 . 0000E+00 
0 . 0000E+00 +,- 0 . OOOOE+OO 
-0 . 1316E-01 +,- 0 . 2000E+00 
-0 . 2518E-03 +,- 0 . 5000E-02 
-0 . 3003E-04 +,- 0 . 1000E-02 
0 . 5494E-03 +,- 0.5000E-02 
-0 . 1634E-03 +,- 0 . 5000E-02 
-0 . 1458E-03 + , - 0 . 5000E-02 
0. OOOOE+OO +,- 0. OOOOE+OO 
0. OOOOE+OO +,- 0. OOOOE+OO 
0. OOOOE+OO +,- 0. OOOOE+OO 
0. OOOOE+OO +,- 0. OOOOE+OO 
0. OOOOE+OO +,- 0. OOOOE+OO 
0. OOOOE+OO +,- 0. OOOOE+OO 
0. OOOOE+OO +,- 0. OOOOE+OO 


ESTIMATION ERROR MEASUREMENT NOISE 


MEAN STDEV 

-0.1833E-01 +,- 0 . 2263E+00 
0 . 1263E-02 +,- 0 . 1888E-01 
-0 . 2547E-03 +,- 0.2949E-01 
-0 . 4934E-05 +,- 0.8413E-04 
-0 . 2701E-04 +,- 0 . 7930E-04 
0 . 1463E-03 +,- 0 . 1915E+00 
-0 . 3046E-03 +,- 0 . 5002E-02 
0 . 2111E-05 +,- 0 . 6532E-03 
0 . 4474E-03 +,- 0.8654E-01 
-0 . 1903E-03 +,- 0 . 1193E-01 
-0 . 7746E-03 +,- 0.7100E-02 
0 . 1017E-03 +,- 0 . 2771E-02 
0 . 2260E-03 +,- 0 . 7848E-02 
-0 . 2969E-13 +,- 0 . 9923E-02 
-0 . 4722E-02 +,- 0.3295E-01 


MEAN STDEV 

0. OOOOE+OO +,- 0. OOOOE+OO 
0. OOOOE+OO +,- 0. OOOOE+OO 
0. OOOOE+OO +,- 0. OOOOE+OO 
0. OOOOE+OO +,- 0. OOOOE+OO 
0. OOOOE+OO +,- 0. OOOOE+OO 
-0 . 1316E-01 +,- 0 . 2000E+00 
-0 . 2518E-03 +,- 0 . 5000E-02 
-0 . 3003E-04 +,- 0 . 1000E-02 
0. 1099E-01 +,- 0 . 1000E+00 
-0 . 6538E-03 +,- 0.2000E-01 
-0 . 5832E-03 +,- 0.2000E-01 
0. OOOOE+OO +,- 0. OOOOE+OO 
0 . 1201E-03 +,- 0 . 1000E-01 
0 . 4505E-03 +,- 0 . 1000E-01 
-0 . 1008E-01 +,- 0 . 5000E-01 


D3.- Summaries of accident solution results (a) problem 1, 

b) problem 2. 
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APPENDIX E 


USER SUBROUTINE 


This appendix describes the subroutine DATA that the user must prepare to 
provide the link between a source of flight data and the SMACK state estimation 
program. Some chores the user may wish to perform in DATA include wild-point win- 
dowing, amplitude scaling, sample compression, and air-data calculations. The user 
may prefer, of course, to do these chores separately and provide a file of time 
histories prepared especially for easy processing by subroutine DATA. Compensation 
for eg offset of accelerometer and air-data measurement systems is optional in DATA 
since position errors ray be specified by including ACC, P-S, and VNE statements in 
the coding list. The first few lines of code for subroutine DATA should appear as 

follows: 

SUBROUTINE DATA (IFIN) 

COMMON/MSM/Z(30,6000) ,Y(30,6000) ,KZN(30,6000) , 

* H(30, 6000), D(30, 6000), S(30,6000) 

C0MM0N/STM/X( 2 1 ,6000) ,W( 9,6000) ,DMW( 54000) 

COMMON/ AUX/XN( 6000) ,YN( 6000) ,ZN( 6000) ,FN( 6000) , 

* KFN ( 6000 ) ,NFN( 6000) 

C0MM0N/SKL/CMK , CMF , CRD , CNL , CMM , CMG 
COMMON / CON /RG AS , TSSL , RH00 , PI , C( 30 ) 

C0MM0N/TIM/DT , NPTS , NI NT , NSMP , NSKP , MPTS , FC 1 

COMMON/MOD/ ALF0,CN AO, WGLD,V ARM 


These arrays and parameters ar<= defined as: 


Z( 30, 6000) 

Flight data array to be filled 

KZN(30,6000) 

Companion array to be filled 

Y(30,6000) 

Output estimates returned with solution 

H(30,6000) 

Output estimates returned with solution 

D( 30, 6000) 

Auxiliary array 

S( 30, 6000) 

Auxiliary array 

X(21 ,6000) 

State-variable estimates returned with solution 

W( 9,6000) 

Forcing-function estimates returned with solution 

DMW( 54000) 

Auxiliary array 

XN(6000) 

Auxiliary array 

YN(6000) 

Auxiliary array 

ZN(6000) 

Auxiliary array 

FN(6000) 

Auxiliary array 
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KFN(6000) Auxiliary array 
NFN(6000) Auxiliary array 

CMK Conversion to meter /sec from knots (0.5144) 

CMF Conversion to meters from feet (0.3048) 

CRD Conversion to radians from degrees (0.01745) 

CNL Conversion to newtons from pounds (4.4482) 

CMM Conversion to meters from nautical miles (1852.) 

CMG Conversion to meter/s 2 from g units (9.807) 

DT Time step for integrating equations of motion, seconds 

NPTS Number of time steps per record (maximum is MPTS) 

NINT Number of integration steps per datum point 

NSMP Number of datum points per record 

NSKP Plotting interval for output 

MPTS Maximum number of time steps per record (6000) 

FC1 Low-pass-filter cutoff frequency 

ALF0 Zero-lift angle of attack, radians 
CNA0 Normal force derivative, /radian 
WGLD Wing loading, newton/m 2 
VARM Magnetic variation, radians 

RGAS Gas constant (286.924 m 2 /s 2 /deg K) 

TSSL Sea-level temperature (288.15 deg K) 

RH00 Air density at sea level (1.225 kg/m^) 

PI Value of Pi (3.14159265 radians) 

C(30) Array of constants set by coding list 

Subroutine DATA is called by SMACK twice, once to fill the appropriate rows of 
the Z and KZN arrays with measurement records, and again after the solution has been 
completed. Entry to DATA is determined by the value of IFIN (0 or 1). Each entry 
must be accompanied by a separate RETURN statement. 

Both the Z and KZN arrays are initialized to zero before SMACK calls DATA the 
first time (IFIN = 0). When a "good" datum point is entered in a row of Z, a 1 
should be entered in the corresponding Location of KZN. In this way channels with 
wild points and/or nonuniform sampling rates can be handled in a mathematically 
correct way (no weight for a bad or missing point). Rows of the flight-data arrays 
Z and KZN are defined in table 5.1. Note that internally, all rotational and linear 

motion variables are in radian ana mks units. The conversion constants CMK, ,CMG 

and the other constants in block /CON/ are provided for user convenience. The 
auxiliary time-history arrays may be used as temporary storage in performing neces- 
sary data-processing chores. 

As indicated in chapter 5 which covers the preparation of the coding list, DT 
is the time step for integration of the state equations and NPTS is the total number 
of points in each state (and output) estimate, not including the initial condi- 
tion. Although these parameters (and NINT, NSMP, and NSKP in block /TIM/) are 
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normally set by the REC statement in the coding list, they may be reset in subrou- 
tine DATA. The user should be aware that each data record (row of Z, KZN) must have 
samples that occur at integer multiples of DT. Nonuniform sample rates are permis- 
sible as long as every interval is divisible oy DT. It should be emphasized that 
fc each "good" datum point entered in a row of Z, a 1 should be entered in the cor- 
responding location of KZN. 

The aircraft-dependent performance parameters ALFO, CNAO, and WGLD in block 
/MOD/ need to be specified in DATA only when an estimate of angle of attack is 
needed but no measurement is available, as in an accident analysis from radar track 
and wind data. For example, for small values the estimate can be approximated by 

AV = ALFO - WGLD * AZ/( CUE*CNA0 ) (El) 

where AV is angle of attack, AZ is vertical acceleration, and CUE is dynamic pres- 
sure. This estimate is provided by SMACK for the radar solution illustrated by the 
first problem of appendix D; in other situations the user should compute the pseudo 
measurement of angle of attack within DATA, when necessary. 

Note that IFIN is zero when DATA is first called. When the solution has been 
completed, IFIN is set to one and DATA is again called. At this time the user may 
wish to perform other operations with the data or the output estimates, which are 
returned in the arrays H and Y. The ith row of the output estimate Y is calculated 

in SMACK as 

Y( I ,N) = SF( I ) *H( I , N ) + VB( I ) , N=1,...,NPTS (E2) 

where SF(I) is the scale factor, VB(I) is the bias error and H(I,N) is the output 
variable, with rows in the order specified by table 5.1. The estimates of aerody- 
namic variables (VT, AV, BV) and specific forces (AX, AY, AZ) stored in array H may 
include positon corrections. The cg-referenced time histories of true airspeed, 
angle of attack, and sideslip angle may be calculated from the air velocities (UA, 
VA, WA) stored in rows 1-3 of array S. The true specific forces (AXB, AYB, AZB) are 
stored in rows 4-6 of array S. 

If it is necessary to compress, filter, or interpolate any records in DATA, the 
user may call subroutine FILT. As indicated in appendix A, this subroutine is used 
in the SMACK starting procedure to provide filtering of each measurement time his- 
tory. The amplitude-frequency characteristic is that of a fourth-order low-pass 
filter with zero phase shift. FILT is based on the backward-filter, forward- 
smoother algorithm used by SMACK, adapted for linear state and measurement models, 
which is discussed in chapter 2 and outlined in table 2.2(b). The filter model and 
its calling statement are described in the following section. 


A Digital Filter 

The low-pass digital filtering subroutine FILT is called several times during 
the SMACK starting procedure. Since the user may have occasion to call FILT from 
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the DATA subroutine, its characteristics and calling statement are discussed here. 
The filter design is based on the solution properties of a fixed-interval smoothing 
problem for a linear, single-input, single-output (SISO), nth-order continuous 
system. It is well known that the solution of this problem yields a filter transfer 
function with 2*n poles equally distributed on a circle centered at the origin of 
the s-plane. The resulting frequency characteristic, therefore, is equivalent to 
two nth-order Butterworth filters in cascade, each with equal bandwidth and equal 
but opposite phase shift at all frequencies. Hence, the filter response is 3*n db 
down at the cutoff frequency, rolls off at 12*n db/octave, and exhibits no phase 
shift. 

The linear continuous problem and its solution are illustrated in figure El for 
a second order system (n=2). Notice that the forcing function is the second 
derivative of the filter output. A finite-difference approximation for this case is 
used with the linear algorithm of table 2.2(b) as the basis for subroutine FILT. It 
has a state model with 
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and a measurement model with 

H = (1 0 0) (E4) 

For an SISO system, the weighting matrices Q and R are scalars: the filter realiza- 

tion determines the ratio to be 


Q/R = (2*f c h) 4 (E5) 

where h is the integration time step (in seconds), and f is the filter cutoff 
frequency (in Hz). The resulting response will be 6 db down at f , and the 
attenuation rate will be 24 db/octave. The extra state in both the continuous and 
discrete state models is added to accommodate an average value of the forcing 
function. It does not affect the second-order nature of the response. The calling 
statement for FILT is 

CALL FILT(FN, KFN, XN, YN, WN, XO, YO, WO, FC, H, NPTS) 
where 

FN Data record to be filtered (destroyed on return) 

KFN Companion array: each "good" point of FN should have a 

corresponding entry of KFN set to one; otherwise the 
entry should be set to zero 
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(a) FIXED-INTERVAL SMOOTHING PROBLEM: 
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STATE MODEL: x 1 =x 2 , x 1 (0) = x 1 q 

x 2 = w + b , x 2 (0) = x 20 
b = 0 , b(0) = b 


MEASUREMENT MODEL: z = x-,+v, FOR (0, T) 

CHOOSE x 1Q , x 2Q , b AND FORCING-FUNCTION w(t) TO MINIMIZE: 

J = (1/2) / [w 2 /q + (z - x-j) 2 /r] dt 
0 


(b) SOLUTION PROPERTIES: 

TRANSFER FUNCTION: 

H(s) = X 1 (s)/Z(s) 

= 0J 4 /(s 4 + OJp) 

WHERE cj c = 27rf c = (q/r) 1/4 
FREQUENCY RESPONSE: 
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Figure El.- Basis for digital-filtering subroutine FILT. (a) Fixed-interval 
smoothing problem, (b) fourth-order frequency characteristic. 
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XN Filtered estimate of FN (array of MPTS) 

YN First derivative of XN (array' of MPTS) 

WN Second derivative of XN (array of MPTS) 

XO Initial condition of XN 

YO Initial condition of YN 

WO Initial condition of WN 

FC Low-pass cutoff frequency, in Hz 

H Integration step, in seconds 

NPTS Number of integration steps (NPTS < MPTS) 

When FILT is used as an interpolating filter, the cutoff frequency FC should be 
set less than or equal to one-half the sampling frequency (H is the data-sample 
interval). If data compression is required after filtering, the final sampling rate 
should be greater than or equal to twice FC. To obtain good results with FILT, the 
original record length NPTS should be greater than 50 points. 
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